General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



TECHNIQUES OF ORBITAL 
DECAY AND LONG-TERM 
EPHEMERIS PREDICTION FOR 
SATELLITES IN EARTH ORBIT 


H 


Prepared for 

NATIONAL AERONAUTICS AND 
SPACE ADMINISTRATION 

GEORGE C. MARSHALL 
SPACE FLIGHT CENTER 

AERO-ASTRODYNAMICS LABORATORY 
HUNTSVILLE, ALABAMA 


/# A A 

0 DEC 1971 'Jt 

(g RECEIVED 
\ 5 , laniMiq 

lS# y 


NOVEMBER 1971 


N 7 2 - 12581 


1 


Unclas 


10063 


BECjy AND LONG-LbI EPHP«f KS ° F ° Bf,ITA '- 

FOR satellite : is r, s n PREU *cri JN 

iJ.’ 1 (Computer Sciences corp tidrl l'« 

_ .141 p *■ '-orp.) Nov. 197] 

5 (NASA CR OK IMA W* r e rr , 

^^Cl£2c G3/2 1 



COMPUTER SCIENCES CORPORATION 






TECHNIQUES OF ORBITAL 
DECAY AND LONG-TERM 
EPHEMERIS PREDICTION FOR 
SATELLITES IN EARTH ORBIT 


Prepared for 

NATIONAL AERONAUTICS AND 
SPACE ADMINISTRATION 


GEORGE C. MARSHALL 
SPACE FLIGHT CENTER 

AERO-ASTRODYNAMICS LABORATORY 
HUNTSVILLE, ALABAMA 


NOVEMBER 1971 


COMPUTER SCIENCES CORPORATION 

8300 South Whitesburg Drive 
Huntsville, Alabama 35802 


Major Offices and Facilities Throughout the World 


TECHNIQUES OF ORBITAL DECAY AND LONG-TERM 
EPHEMERIS PREDICTION FOR SATELLITES IN EARTH ORBIT 


November 1971 
by 

B. F. Barry 
R. S. Pimm 

C. K. Rowe 

Prepared For 

National Aeronautics and Space Administration 
George C. Marshall Space Flight Center 
Aero-Astrodynamics Laboratory 
Contract NAS8-26113 


APPROVED BY: ^ ^ 

W. B. Twichell, Jr. 
Director 

Aerospace Systems Center 



COMPUTER SCIENCES CORPORATION 
8300 South Whitesburg Drive 
Huntsville, Alabama 35802 


FOREWORD 


This report, which supersedes the interim report dated May 1971, presents the 
results of work performed by Computer Sciences Corporation's Aerospace Systems 
Center while under contract to the Aero-Astrodynamics Laboratory of the George C. 
Marshall Space Flight Center, Contract NASS -2G 113. 

The authors are grateful to Messrs. L. D. Mullins and B. S. Perrine (MSFC- 
S&E-AERO-MMD) for their technical assistance and to Messrs. W. J. Elkins and 
M. M. Hansing (CSC) for their programming support. 


SUMMARY 


The methods of both special and general perturbation theory are employed in 
solving the equations of motion for a satellite subjected to the perturbational effects 
of earth oblateness and atmospheric drag. In the special perturbation method, 

Cowell and variation-of -parameters formulations of the motion equations are imple- 
mented and numerically integrated by means of a MARVES (Marshall Vehicle Engi- 
neering Simulation System) computer program. Variations in the orbital elements 
due to drag are computed using the 1970 Jacchia atmospheric density model, which 
includes the effects of semiannual variations, diurnal bulge, solar activity, and geo- 
magnetic activity. In the general perturbation method, two-variable asymptotic 
series and the automated manipulation capabilities of FORMA C (Formula Manipulation 
Compiler) are used to obtain analytical solutions to the variation-of-parameters 
equations. Solutions are obtained when considering the effect of oblateness only 
(J and J, ) and the combined effects of oblateness and drag. These solutions are 

u t) 

then numerically evaluated by means of a FORTRAN program in which an updating 
scheme is used to maintain accurate epoch values of the elements. The atmospheric 
density function is approximated by a Fourier series in true anomaly, and the 1970 
Jacchia model is used to periodically update the Fourier coefficients. The accuracy 
of both methods is demonstrated by comparing computed orbital elements to actual 
elements (or elements computed by standard MSFC programs) over time spans of up 
to 8 days for the special perturbation method and up to 356 days for the general 
perturbation method. 
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NOMENCLATURE 


Mathematical Symbols 

a semi major axis 

a. b f c a set of constants arising in the solution of the ordinary 

differential equations having t_ as the independent variable 
(explicitly defined in Appendix F) 

a.,b. Fourier coefficients appearing in Fourier series 

approximation to atmospheric density function 

A a constant arising in the solution of the ordinary differential 

equations having Tas the independent variable (explicitly 
defined in Appendix E) 

(A/m) satellite arca/inrss ratio 

b angle measured normal to orbital plane in direction of normal 

perturbative acceleration 

-1 /2 

B orbital element defined as a 

C r C ... a set of constants arising in the solution of the ordinary 
differential equations having t_as the independent variable 
(explicitly defined in Appendix E) 

integration ’’constant" associated with asymptotic series 
solution development (see Paragraph 4,3.1) 

aerodynamic drag coefficient 

drag force magnitude per unit mass 

a constant arising in the solution of the ordinary differential 
equations having t_ as the independent variable (explicitly 
defined in Appendix F) 

eccentricity 

specific angular momentum 


C(t) 


D 


D 


D„ 
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NOMENCLATURE (Continued) 


inclination relative to earth equatorial plane 

coefficients of second, third and fourth harmonics, respectively, 
of earth gravitational potential 

constant arising in the formulation of the differential equations 
of motion for tangential atmospheric drag (explicitly defined in 
Paragraph 2. 5. 2) 

mean anomaly 

, r- J / 3 1/2 

mean motion defined as (M/a ) 

2 

semilatus parameter defined as a(l-e ) 

geocentric radius vector 

perturbative acceleration vector 

equatorial radius of earth 

perturbative gravitational potential function 

time 

fast time variable defined as t (l+ct £) 

£ 

slow time variable defined as et 
argument of latitude (P + u') 
inertial velocity vector 
relative velocity vector 

a constant arising in the solution of the ordinary differential 
equations having T as the independent variable (explicitly 
defined in Appendix E) 



NOMENCLATURE (Continued) 


ol 8 v 6 
V i’ 7 f i 


a 

a 

f 

/?' 

6 

c 

V 

\ 


V* 


4 

p 

0 

a; 

w 


n 


a set of constants used to represent different linear combinations 
of the Fourier coefficients a. and (explicitly defined in Appendix F) 

right ascension of satellite subpoint 

flight path angle (positive above local horizontal) 

angle between local latitude and orbital planes 

declination of satellite subpoint 

small perturbative parameter defined as (:</2)<J 2 

transformation parameter defined as e since 

angle between local longitude and orbital planes 

constants appearing in the £ and 7£ solutions for oblateness/ 
drag (explicitly defined in Appendix F) 

earth gravitational constant 

true anomaly 

transformation parameter defined as e cost e 
atmospheric density 

angle between radius and velocity vectors 
argument of perigee 

magnitude of earth rotational velocity vector 
right ascension of ascending node 


NOTE: The subscript "0" is used to denote the epoch (or reference) value of an 
element or element function (i.e. , e^, e^\ e^D). 
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NOMENCLATURE (Continued) 


Element Typos and Variations 

long-periodic variation - A variation periodic with respect to C£or 
multiples of (£; for example, sino,'. 

mean orbital elements - The osculating elements with the short periodic 
variations removed. 

osculating orbital elements - The instantaneous elements defining the 
continually changing elliptical orbit. 

secular variation - A steady nonose illatory variation from the epoch 
value, i.e., a variation directly proportional to the independent 
variable; for example, Ct. 

short-periodic variatio n - A variation periodic with respect to linear 
combinations of V and UJ; for example, cos (t'+co). 


x 



SECTION 1 - INTRODUCTION 


The purpose of this study is to develop techniques of orbital decay and long- 
term ephemeris prediction for satellites in elliptical earth-orbit. These techniques 
are to be accurate and flexible, and are to lead themselves to rapid computation. 

In order to meet current needs, emphasis is to be placed on the development of 
ephemeris prediction techniques for low-eccentricity, near-earth orbits when con- 
sidering the perturbational effects of earth oblateness and atmospheric drag. 

Classically, two general methods of attack are available for solving this 
problem. These methods are known as the methods of special and general pertur- 
bations, respectively. In both methods, the equations of motion may be formulated 
either as three second-order differential equations (for the perturbative accelerations) 
or as six first-order differential equations (for some set of fundamental orbital ele- 
ments). The two methods differ in that special perturbation formulations (such as 
Cowell's, Encke's, varia*ion-of -parameters, etc.) employ various numerical inte- 
gration procedures (such as Runge-Kutta, Fehlberg, Shanks, etc.) to obtain the 
solution, while general perturbation techniques (such as variation-of -parameters, 
variation-of-coordinates, etc.) generally employ series expansions (such as Taylor's, 
multivariable asymptotic, etc. ) combined with analytical integration to achieve the 
desired solution. In choosing one method or the other, one must keep in mind both 
the nature of the orbit under consideration and the nature of the solutions desired. 

The main advantages of the special perturbation method lie in simplicity of 
formulation, applicability to any type of orbit in any perturbing force field, and com- 
pactness of storage requirements for program solution. This method is ideally 
suited for calculating orbits of limited duration. The main disadvantages inherent 
in this method are the inducement of errors (truncation and round-off) due to the 
numerical nature of the process, the resulting lack of application to orbits of long 
duration, and the extensive computation time required for solution. 
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The primary advantages of the general perturbation method lie in its applica- 
bility to orbits of long duration, its relatively rapid computer solution time, and its 
ability to provide a clearer geometric conception of the effects of the various pertur- 
bations. On the other hand, in applying this method one is faced with much analytical 
labor in formulating the equations to include various perturbations and in obtaining 
the solutions to these equations. 

To achieve extended applicability in attacking the problem at hand, it was 
decided to employ formulations of both methods. In the special perturbation method 
both the Cowell and the variation -of -parameters formulations are employed, while 
the general perturbation method consists of the variation-of -parameters formulation 
using two -variable asymptotic series expansions. To alleviate the analytical labor 
required, the automated manipulation capabilities of the FORMA C (Formula Manipu- 
lation Compiler) language are utilized. 
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SECTION 2 - DERIVATION OF THE VARIATION-OF-PARAMETERS 
DIFFERENTIAL EQUATIONS OF MOTION (LAGRANGE’S 
PLANETARY EQUATIONS) 


The purpose of this section is to derive, by the method of perturbative differ- 
entiation, the differential equations of motion for a selected set of orbital elements 
(or parameters) when considering the perturbational effects of earth oblateness and 
tangential atmospheric drag. Since perturbative forces are additive, the differential 
equations for each perturbational effect can be formulated separately. This set of 
differential equations will then be solved numerically by the methods of special pertur- 
bation theory and analytically by the methods of general perturbation theory. 

2. 1 SELECTED ORBITAL ELEMENT SET 

The orbital element set selected for consideration is 


where B 
a 
e 
i 

Q 

CO 

M 

V 

t 


(B, e, i, O, co, M(or v); t) 

_l/2 

a (defined for mathematical simplification) 

semimajor axis 

eccentricity 

inclination relative to earth equatorial plane 
right ascension of ascending node 
argument of perigee 
mean anomaly 
true anomaly 

time (independent variable) 


Although only one anomaly angle is needed in the element set, it is advantageous 
to consider both M and v. The differential equation for M is more amenable to 
asymptotic series solution; on the other hand, it is mathematically easier to derive the 
differential equations for all elements in terms of u. Consequently, the differential 
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equation for M is derived and solved; V is then obtained by a Fourier-Bessel expansion 
involving M and e. 

2. 2 PERTURBATIVE DIFFERENTIATION 

In the theory of perturbative differentiation, the variation (time -derivative) of 
any element Ms considered as the sum of two parts; i. e. , 

— = f + f x 
dt 

where f_ (f-dot) is the Keplerian variation that remains if all disturbing forces are 
suddenly removed and f^ (f-grave) is the perturbative variation caused by the disturbing 
forces. There are three types of variations which arise in the theory; namely, 

r — - f, where f = 0 

Type 2 : df \ • 

— - f, where f - 0 

Type 3 : df • \ , 

— - f + f, where both parts exist 


Since the velocity associated with the osculating orbit at the point of tangency is 

dr 

the same as the actual velocity of the perturbed satellite, the components of — in an 
inertial coordinate system are of the first type. Variations of the second type arise 
for elements that would be constant in Keplerian motion, such as a, e, i_, Q_, and ox 
Elements referred to a perturbed reference direction, such as and y, are of the 
third type. 

It follows, then, that the basic differential equations of motion for the selected 
elements are 


dB 

dt 


b' 


( 2 - 1 ) 


de 

dt 


( 2 - 2 ) 
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di 

dt 


1 


(2-3) 


dO 

dt 


O' 


(2-4) 


dec 

dt 


= or 


(2-5) 


dM 

dt 


M + M V = n + M x 


( 2 - 6 ) 


The next step is to obtain the perturbative variations indicated above. 

Two techniques are commonly used to obtain the perturbative variation f^ of an 
element f. The first technique consists of developing the total variation of the element 
and then removing the Keplerian part; i. e. , 


f '= 


df 

dt 


f 


The second technique consists of using perturbative differentiation, which involves 
taking the grave -derivative of a given expression in which only the variations due to 
the disturbing forces are considered. The second technique is used here to obtain the 
perturbative variations of the elements. (For a further discussion of perturbative 
differentiation, see References 1 and 2; particularly p. 21 of Reference 1.) 

2. 3 PERTURBATIVE VARIATION EQUATIONS 

It is necessary to obtain the perturbative variations f£, e\ i\ O', oj' and m' in 

±.\ 

terms of the orbital elements and the perturbative acceleration vector r, resolved as 
follows (Reference 2, p. 284): 

r'= r' U + r£'v + rOw 

where U = unit vector in direction of increasing r (radial) 

V = unit vector perpendicular to r in orbital plane (transverse) 

W = unit vector perpendicular to orbital plane (orthogonal) 


2-3 


As will be seen in the next section, the perturbative acceleration components r\ r 
and rb'can also be obtained in terms of the orbital elements via the disturbing function 
R. 

Although derivation of the perturbative variation equations by the method of per- 
turbative differentiation is straightforward, it is mathematically tedious; consequently, 
the procedure is presented in Appendix A. The results are (also see Reference 1, 
p. 22, and Reference 2, pp. 247 and 284): 


b' 

e' 


Q' 



2 k' 

r b cos u 

v/W 


^ - 
r b sin u 

ftp sin i 


(2-7) 


( 2 - 8 ) 

(2-9) 


( 2 - 10 ) 


= 


M' = 


-O' cos i - ~7 = — cos + i\ sin v 

ft? e \r ) ftp e \r / 


-(!-/) 


2\l/2 


U)' + o' cos i + 


2r_r' 

ftp 


( 2 - 11 ) 

( 2 - 12 ) 


2.4 PERTURBATIVE ACCELERATION COMPONENTS 


2. 4. 1 Earth Oblateness 


•> 


The perturbative acceleration vector r due to an axially symmetric oblate earth 
can be written as the gradient of the perturbative potential function R (per unit mass), 
which becomes, in spherical coordinates, 

,o 


* \ = dR -r f 1 d&j 1 5R j- 
r dr 1 r cos 6 /fa' ■' r d6 



where (see Figure 2-1) 


i - unit vector in direction of increasing r 

J - unit vector in direction of increasing a' 

k = unit vector in direction of increasing 6 

and (Reference 3, p. 49) 



(NOTE: Jg and are negative numbers. ) 


As previously mentioned, the general expression for the perturbative acceleration 
vector can be written as 


r =» r v U + rp V + rfi'w 


where U 
V 
W 


unit vector in direction of increasing r (radial) 

unit vector perpendicular to r in orbital plane (transverse) 

unit vector perpendicular to orbital plane (orthogonal) 
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The transformation between the (i, j, k) system and the (U, V, W) system is 
obtained via a right-hand rotation about the i-axis through an angle , as seen 
from Figure 2-1. Thus 



resulting in the scalar equations 


= dJL 

dr 

- sin j3~ 3R 
r 36 

= cos p ' 3R 

r "5F 

Performing the indicated differentiation on Equation (2-13) yields 
r' = - | #1 J 2 r e 2 | ^ (1-3 sin 2 5 ) -2m J 3 r e 3 ^ j (3-5 sin 2 6) sin 6 + ^ m J 4 r e 4 ^ ^ J (35 sin 4 6-30 sin 2 6 + 3) (2-14) 

n = sin p . 3 m J 2 T e^\ j s * n ^ cos 3 + 3 3 r e^‘g s * n ^ 6 )cos 5 — ^ m J 4 r e ^| g j (7 sin 2 6—3) sin 6 cos sj (2—15) 

rb' = cos(3'j^— 3 m J2 r e ^j sin 6 cos 5 + ^ ^ J 3 r e 3 ^ j(l— 5 sin 2 6) cos 6 — "^M J4 r e 4 ^ — ^(7 sin 2 6 3 )|sin 6 co* ( 2 - 16 ) 


f ’ 
r v' 
rb’ 


It is now necessary to express and 6^ in terms of the orbital elements. Referring 
to Figure 2-1, from the spherical triangle ABC 

• sin 6 * sin i sin u (2-17) 


cos 6 


V 1— sin^ i sin^ u 


(2-18) 
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Also, from the same triangle 


or 


cos u - ctn X ctn i 


But 


tan X 


ctn i 
cos u 


1 + tan^X = — - — 
cos^ X 


thus 1 1 cos u 

n/i uan^ X >/ 1 + / c - fo - 1 \J ctn^ i + cos^ u 
\ cos u f 

Since $ 9(T - \ (i. e. , latitude and longitude lines are perpendicular), 

then 


sin /3 ' = sin (90°— X) = cos X 


or 


and 


sin/3' - 
cos |3 ' - 


cos u 


V ctn^i + cos^ u 

ctn i 

V ctn^ i + cos^ u 


However, 

■J ctn^ i + cos^ u = y/ cos^ i + cos^ u sin^ i = -r^-r J 1— sin^ i sin^ u 

sm l smi v * 

hence 


sin /3 ' ~ 


cos/3' = 


cos u sm i 


\/ 1— sin^ i sin ^ u 


cos i 


\J 1— sin^ i sin^ u 


(2 


(2 


Substituting Equations (2-17) through (2-20) into Equations (2-14) through (2-1S) 
yields, after simplification, the desired results. 


- 19 ) 

- 20 ) 
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TV 


rb' 


= J 2 r e 2 (; 

^Vl— 3 pin 2 i sin 2 u) — 2/u Jg r e 3 (3—5 sin 2 i sin 2 u) sin i sin u 


+ 1 (,J 4 r . 4 (. 

' \r / 

(35 sin^ i sin ^ u —30 sin 2 i sin 2 u + 3) 

(2-21) 

= - ! M .1 2 r e 2 ( 

-^jsin 2 i sin 2 u + -| ju Jg r p 3 ^ (1—5 sin 2 i sin 2 u) cos u sin i 

sin 2 i sin 2 u (7 sin 2 i sin 2 u —3 ) 


“I ^ J 4 r e 4 ( 

(2-22) 

* ~ 3 " J 2 r e 2 ( r l 

j sin i cos i sin u + ^ n Jg r g 3 ^ j (1—5 sin 2 i sin 2 u) cos i 

(2-23) 

" J 4 r e 4 ( 

~ ^ sin 2 i sin u (7 sin 2 i sin 2 u — 3) 



After converting to units of earth-radii and performing trigonometric-identity 
manipulations, it can be shown that Equations (2-21) through (2-23) agree with 
Reference 2, p.288, and Reference 4, p. 193. 

2.4.2 Tangential Atmospheric Drag 

The perturbative acceleration vector f'due to a tangential atmospheric drag 
force can be written as 

r ~ - DT + ON + OW - -DT 


where (see Figure 2-2) 

T = unit vector along orbit tangent in direction of motion 
(tangential) 

N - unit vector perpendicular to orbit tangent (normal) 

W = unit vector perpendicular to orbital plane (orthogonal) 



C D p V R 2 
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2 -: 



The relative velocity v can be approximated in terms of the inertial velocity 

R 


vby (Reference 5, p. 165) 


cos i \ 

1 — — ) 


Since the inertial velocity of a satellite in an elliptical orbit is given by (Refer- 
ence 6, p. 80) 

2 * - (1 + e 2 + 2e cos v ) = — (l + e 2 + 2e cos v) 
v P (1— e 2 ) 


the drag force magnitude per unit mass can be written as 


D = i(^) C D P (I~2) (1 + e2 + 2e 005 * ( 1_ 


(2-24) 


As discussed in the previous section, the general expression for the perturba- 
tive acceleration vector is 


r = r ' U + rv' V + rb ’ W 

From Figure 2-2, it can be seen that the transformation matrix relating the 
(U, V, W) system to the (T, N, W) system is obtained via a right-hand rotation about 
the W-axis through an angle (180° + 0), i. e. , 

cos ( 180 ° + 0) sin ( 180 ° + 0) o\ 

0 I = 

1 / 


[T] = ( -sin (180° + 0) cos (180° + 0) 



Thus, 



resulting in the scalar equations 


r 

TV 

rb 


D cos 0 
— D sin 0 

o' 


(2-25) 

(2-26) 

(2-27) 


2-11 


The angle 0 is related to the orbital elements by (Reference 6, p. 83) 

1 + e cos v 


sin 0 = 


cos <p = 


(1 + e 2 + 2e cos v) */ 2 


— e sin v 


(l + e 2 + 2e cos v) ^ 2 


( 2 - 28 ) 

(2-29) 


Substituting Equation (2-24) and Equations (2-28) and (2-29) into Equations (2-25) 
through (2-27) yields the desired results 


1 /A \ ^ e sin v 

r ' =4 (s) c d'’ 


(1—e 2 ) 


(1 + e 2 + 2e cos v) 1/2 (l- 




juB^ (1 + e cos v) 
( 1—e 2 ) 


(l + e 2 + 2ecos..)l/2 2 


(2-30) 


(2-31) 


rb' = 0 


(2-32) 


2. 5 DIFFERENTIAL EQUATIONS IN FINAL FORM 
2.5.1 Earth Oblateness 

Expressing the earth oblateness differential equations in final form requires 
substituting Equations (2-21) through (2-23) into Equations (2-7) through (2-12), 
simplifying, and then substituting the corresponding results into Equations (2-1) 
through (2-6). To illustrate this procedure, the final form of the differential equation 
for the element i__ will be derived. The equations for all other elements can be 
obtained in a similar manner. 

Substituting Equation (2-23) into Equation (2-9) yields 


i ’ 


cos 


n/mp 


f [-3pJ2 r e 2 (^) sini 


i cos i sin u + 


^ M «J 4 r e ^ ^ j sin 2 i sin u (7 sin 2 i sin 2 



(l — 5 sin 2 i sin 2 u) 


cos i 
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Substituting 


r = 


a(l— e 2 ) P 


1 + e cos v l + e cos v 


and performing trigonometric -identity manipulations yields 

3. v 4rj 2 r e 2 ( 1 + e c°s^) 3 . n . . n 
i = — sin 2 l sin 2 u 


4 P 


7/2 


3\/m’ J 3 r J* (! + e cos y ) 4 o o 

+ j— - cos i cos u (1— 5 sin' 1 i sin u) 

2p 9/2 

5v^TJ 4 r e 4 (l + ecosp) 5 0 0 

^ YYj 2 sin 2 i sin 2 u (7 sin" 2 i sin" 2 u — 3) 


This equation can be rewritten as 

(| ** 2 ) V^r e 2 (l+ecosv X J' 




2p 

Jm 


7/2 


sin 2 i sin 2 u 


(| J 2)(lo) v7r r «“ 3 (1 + e cos " )4 , 2 

+ — ZJA—ttl — ___ cos j cos u (I—5 sin^ i sin z u 

p y /2 

5 (I J2 )(j^) r e 4 (1 + e cos p > 5 


Defining 


12 p 11 / 2 

3 

e = "2 J 2 (a small parameter) 


sin 2 i sin 2 u (7 sin 2 i sin 
< 2 - 


This equation becomes 


Vm r_ 2 (1 + e cos v)' 


1 = — e 


2 P 


7/2 


sin 2 i sin 2 u 


( ^2 ) V^ r e 3 + e cos y ) 4 
+ 6 79/2 


cos 


i cos u (1—5 sin 2 i sin 2 u) 


5 (j^) ^ r e 4 (1 + ecosy) 5 

— e ~~ sin 2 i sin 2 u (7 sin 2 i sin 2 


i2 P n/2 


* u- 3) 
33) 


u-3) 
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Finally, substituting this equation into Equation (2-3). noting that 

I = 1 = B 2 

p a(l-e 2 ) (1 — e 2 ) 

yields the desired result 


di 

dt = ~ e 


\/v~ B 7 r g 2 (1 + e cos v)^ 
2(l-e 2 ) 7 / 2 


sin 2 i sin 2 u 


r H B® r e ^ ( 1 + e cos r ) 4 


+ e 


(1 e ) 


2,9/2 


9 9 

cos i cos u (1~5 sin- 6 i sin u) 


— € 



y//T B 11 r e 4 (1 + e cos v) 5 
12 (1— e 2 ) 11 , 2 


sin 2 i sin 2 u (7 sin 2 i sin 2 u — 3) 


(NOTE: An optional formulation of this equation would be one in which the higher 

earth harmonics (<T and J are treated as ’'higher-order" perturbing terms; i.e. , 

2 

as £_ terms. As will be shown in the solution procedure of Section 4, however, 
treating these harmonics as £_-order terms yields mean orbital elements which 
include long-periodic as well as secular variations. ) 

As previously mentioned, the equations for the other elements can be obtained 
in a similar manner. The complete set of differential equations when considering 
earth oblateness is presented on the following pages. 
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Differential Equations of Motion When Considering Earth Oblateness 


to 

I 


dB 

dt 


\/m~ B 3 r 2 ( l + e cos v)* 


(1-e 


1 + e cos vr f . 9 0 o t 

2 ^ [ e sin u < 1 ~ 3 s,n > sin^u) + (1+e cos i>) sin 2 i sin 2 u 


(2-34) 


+ 4- e 


u 1 


r e ®(l+ecos»*)® 


(,-,. 2 , 11/2 

I J, 


[ e sin p sin 1 sin u (3-5 sin 2 i sin 2 u) -i- (1+e cos v) ( 1-5 sin 2 i sin 2 u) cos u sin ij 


y/J*- B 42 r^ll+ecosi') 


- p 2v13/2 


(1— e^) 


~ [ e SU1 v (35 sin4 * sir * 4 u_30 s *n 2 « sin 2 u-3) - 2 (1 +e cos u) sin 2 i sin 2 u (7 sin 2 i sin 2 u-3) J 


de 

dt 


_ _ s/vB 1 rj 2 (1+e cos v) 


(1-e 2 ) 7 / 2 I- 8 ™ V (1 +C TOS U) (1 ~ 3 Si " 2 * si ” 2 u) + sin2 i si n 2u (2 cos v + e cos 2 v + e) j 


(2-35) 


— — e 


+ -5_ e 

12 


V ^ rb9 (^) r e 3 < 1+ecos 0 4 


( l_ e 2)9/2 
*2^ 

( 1 ^ 2 , 11/2 


[sin i sin u sin v (1+e cos v ) (3-5 sin 2 i sin 2 u) — J-(l-5 sin 2 i sin 2 u) (2 cos v+e cos 2 v+e) cos u sin ij 


y/V r e 4 (l + e cos 


[sin v (1+e cos j) (35 sin 4 i sin 4 u - 30 sin 2 i sin 2 u + 3) - 2 sin 2 i sin 2 u (7 sin 2 i sin 2 u — 3) 

(2 cos v + e cos 2 v + e)J 


ayf 
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(2-36) 


di _ 
dt “ 


— e 


+ e 


x/TTB 7 r 2 (1 +e cos p) 3 

— — sin 2i sin 2u 

2(l-e 2 ) 7 / 2 

r e 3 (1+e cos y) 4 


(1-e 2 ) 9 / 2 


cos i cos u (1 — 5 sin 2 i sin 2 u) 


_5_ 

12 


e 


V^Tb 11 r e 4 (1+e cos 

a-e 2 )H/ 2 


sin 2i sin 2u (7 sin 2 i sin 2 u — 3) 


d« 

dt 


= _o V^B 7 r e 2 (l+e COSI >) 3 9 

e ^z^772 sin ucosi 

\/m"B 9 r e 3 (1 +e cos vf 


(2-37) 


(1— e 2)9/2 
• 1 1 / Jy| 


ctn i sin u (1—5 sin 2 i sin 2 u) 


^ y/fj. r e 4 (1+e cos *>) 3 


(l-e 2 )H/ 2 


cos i sin 2 u (7 sin 2 i sin 2 u — 3) 
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dco _ ^ \/^B 7 r e 2 (l+ecosf ) 3 


, (l+e cos t')' J f 2 2 

^ g2)7/2 [ 2esin u cos^ i + cos i* (l+e cos v) (1—3 sin 2 i sin 2 u) — sin v sin 2 i sin 2u (2+e cos t , )| 


>/? r9 ^) r e 3 ( 1+e cos v ) 4 r 


_ p 2v9/2 


e ( 1 — e" 6 ) 




J^(l 5 sin 2 i sin 2 u) (— e ctn i cos i sin u + sin v cos u sin i j 2+e cos v ! ) + -4. (1 +e cos v) 

I * 3 

(3—5 sin 2 i sin 2 u) cos v sin i sin uj 


r* (1 +e cos 


- - (7 sin 2 i sin 2 u— 3) (e cos 2 i sin 2 u— 4-sin 2 i sin 2 u sin v |2+e cos v j ) — -L(l+ e cos i>\ 

e ( 1 — e*) 11 '* L 2 ( | 4 ' 

(35 sin 4 i sin 4 u — 30 sin 2 i sin 2 u + 3) cos I'J 


dM _ v/^B 7 r 2 (l+e cos v ) 3 


V M r e (l+e cos V) 9 . . 9 . . _ I i „ 

e ^_ e 2 j 3 3sm Isin u ) (~ 2 e+cos v j l+e cos v j ) — sin v sin 2 i sin 2 u ( 2 +e cos r) 

V^B 9 (ii) ^(Hecos,! 4 

' " 2 ' r „ , , 


e ( 1 — e 2 ) 4 


x/m” B 44 ( ™) r e 4 (l+e cos v)® 


e (i_ e 2)5 


£ sin i sin u (3 5 sin 2 i sin 2 u) (— 2 e+cos v j 1 +e cos v J ) +-3-sin v cos u sin i ( 2 +e cos i>) (1 —5 sin 2 i sin 2 u)j 

£(35 sin 4 i sin 4 u 30 sin 2 i sin 2 u + 3) ( — 2e + cos v jl+e cos v j ) + 2 sin 2 i sin 2 u (7 sin 2 i sin 2 u — 3) 

sin v ( 2 +e cos v)J 



2,5.2 Tangential Atmospheric Drag 

Expressing the drag differential equations in final form requires substituting 
Equations (2-30) through (2-32) into Equations (2-7) through (2-12), simplifying, and 
then substituting the corresponding results into Equations (2-1) through (2-6). To 
illustrate this procedure, the final form of the differential equation for the element 
e will be derived. The equations for all other elements can be obtained in a 
similar manner. 


Substituting Equations (2-25) and (2-26) into Equation (2-8) yields 


= \fi-P [ fD C ° S * (^ sin ~ rD sin + *) cos v * « j ] 


Substituting 


r = 


1+e cos v 


JL * 1+e cos v 
r 


the equation becomes 

D r p cos 0 


Vfp t 1+e 


cos V 


p sin 0 pesim 

(1 +e cos v) sin v — - — (2+e cos v) cosy — 


1+e cos v 


1+e cos v 


D\Zp~ 

y^l+e cos v) 


[cos 0 sin v (1 +e cos v) — sin 0 cos v (2+e cos v) — e sin v] 


Substituting Equations (2-28) and (2-29) yields 


Qy'p 


' = — v * ~ n i in [— e sin 2 v (1+e cos v) — cos v (1+e cos v) (2+e cos v) -e (1+e cos y)] 

i/iMl + e cos y) (1+e + 2e cos y) ' 


2Dv/p^ 

= — — 9 j-pt (e + cosy) 

+ 2e cos y) 1 ' 2 

Substituting Equation (2-24) and noting that 


p«a(l— e 2 ) -i 1 . * 5 


results in 


2y/l— e 2 (e+cos y) 

B v ' r iT(l+e 2 + 2ecosy)^ 2 


[*(=) 


C D P {1 _ e 2j ( 1 +e 2 + 2e cos v) 


](- 


(A.) 

t Cq py/jTE (e+cos u) (1+e 2 + 2e cos y)*/ 2 j 

f , _ w e cos i \ 

\m ) 

' (l-e 2 )l/ 2 1 

l n 1 
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Since 


n=V>TB 3 

the equation can be written as 

< m lA\ ^ (e+<»8 v) (1 +e 2 + 2e cos v)*^ 2 

W „(l- e 2)l /2 




In terms of the small parameter £ as defined by Equation (2-33), this becomes 

/ A \ / 1 \ 2 C D p H B 4 (e+cos v) (1+e 2 + 2e cos r) 1 / 2 / o; e cos i 
\ m M«V 3n(l-e 2 ) 1/2 V n 

Finally, substituting this equation into Equation (2-2), and defining 

; \ 2 


K* = 




yields the desired result 


_de = K* p (e+cos p) (1+e 2 + 2e cos m)*/ 2 
dt n (l- e 2)l/2 


(2-40) 


As previously mentioned, the equations for the other elements are obtained in a 
similar manner. The complete set of differential equations when considering 
tangential atmospheric drag is presented on the following page. 
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Differential Equations of Motion When 
Considering Tangential Atmospheric Drag 


dB _ € K* p ( 1 +e 2 + 2e ccs v ) 3 / 2 

dt 2n(l-e 2 ) 3/2 


(2-41) 


de = 
dt 


— e K* P (e+cos y) (1 +e 2 + 2e cos v 


n (1— e 2 )^^ 2 


(2-42) 



(since r b 0) 


(2-43) 


,$£L = 0 (since rb' = 0) 


(2-44) 


= — e K* p B^ sin v (l+e^ K 2 eco S ^/2 
ne(l— e 2 )^/ 2 


(2-45) 


dM 

dt 


e K *P B 4 sin y(l+e 2 + 2e cosy) 1 / 2 j \ + e ’ 
n \e 1+ecost/j 


(2-46) 
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SECTION 3 - SPECIAL PERTURBATION METHOD OF SOLUTION 


In Section 2, the differential equations of motion for a selected set of orbital 
elements were derived when considering the perturb at ional effects of earth oblateness 
and tangential atmospheric drag. The purpose of this section is to discuss the pro- 
cedure by which these equations, along with a Cowell formulation of the motion 
equations, are numerically solved. Included is a synopsis of the MARVES computer 
program which has been implemented to perform this so-called special perturbation 
method of solution. 

3. 1 GENERAL 

The term "special perturbations" refers to a technique for the prediction of an 
orbit by numerical integration, so as to include the effects of various perturbative 
forces that cause the trajectory to deviate from some reference orbit (Reference 2, 
pp. 227-228). The basic procedure is the generation of the next step or increment 
of the state variables representing the orbiting body when having a complete knowledge 
of the preceding variables (Reference 7, pp. 220-221). Specifically, one begins with 
some epoch state and integrates, numerically, a set of three second-order or six 
first-order differential equations of motion. 

The variation-of-parameters formulation involves the integration of six first- 
order equations (often referred to as the Lagrange planetary equations) which are 
functions of the selected orbital elements. As is evident, in the literature (Reference 
2, p. 243; Reference 8, pp. 235-230), there is no "best" set of fundamental elements 
to employ, and the choice is dictated by the application in mind. In the Cowell formu- 
lation, three second -order motion equations for the perturbative rectangular accelera- 
tions are integrated to obtain the current state variables (position and velocity). 

3.2 VARIATION-OF-PARAMETERS FORMULATION 

In the variation-of-parameters formulation, six first-order element rate equa- 
tions are numerically integrated; these equations reflect perturbations due to earth 
oblateness (second, third, and fourth harmonics) and atmospheric drag (using a 1970 
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Jacchia atmospheric density model). The oblateness equations are identically those 
presented in Paragraph 2.5. 1, whereas the drag equations differ from those presented 
in Paragraph 2. 5.2 in that p three-dimensional rather than tangential drag force is 
considered. (NOTE: To readily obtain analytical solutions to the equations, it is 
necessary to assume a tangential atmospheric drag; however, this assumption is not 
required when numerically integrating the equations. ) 

3.3 COWELL FORMULATION 

In the Cowell formulation, the equations of motion are expressed in rectangular 
form and integrated twice to obtain the velocity and position. These equations have 
the standard form: 


.2 

= 5 ? £ 

dt 2 


<x-*y, z) 


where x represents the central force term, and x, the perturbative term, represents 
the accumulated effects of all perturbations acting. The perturbations included in this 
formulation also consist of earth oblateness (second, third, and fourth harmonics) and 
atmospheric drag (using a 1970 Jacchia atmospheric density model). 


3.4 SYNOPSIS OF MARVE S COMPUTER PROGRAM FOR NUMERICALLY INTE- 
GRATING THE EQUATIONS OF MOTION 

A MARVES/FORTRAN double-precision special perturbations program has been 
developed for the UNIVAC 1108 and is currently available through the MSEC Computa- 
tion Laboratory. This program provides, on user option, either the Cowell or variation- 
of-paraineters formulations. 


The program is modular in design, with FORTRAN subprograms selectively 
linked and controlled by two MARVES driver programs. This configuration allows 
user selection from a library of simulation routines and high precision numerical 
integration schemes currently operational and available to MARVES users (Reference 
9). These integration schemes include a variety of single and multistep methods with 
provisions for optimum step-size prediction based on the resultant truncation error. 
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When selecting the "best" method for solving the set of differential equations for a 
particular orbit, many factors must be considered, such as: required accuracy, 
number of integration interrupts, frequency of computer printout, and integration 
step-size limitations. Reference 9 contains a thorough discussion of the methods 
currently available in MARVES, along with some generalizations that can be made 
about method selection. 

Many other desired features are incorporated into the program, such as critical 
time events, the nearest Besselian year coordinate transformation, and the 1970 Jacchia 
atmospheric density model. Also included is a solar-ephemeris computation routine 
that eliminates the need for read/interpolation of Jet Propulsion Laboratory (JPL) 
ephemeris tapes. 

A complete description of this MARVES program (referred to as the SPERTB 
program) is given in Reference 1C. 
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SECTION 4 - GENERAL PERTURBATION METHOD OF SOLUTION 


In Section 2, the differential equations of motion for a selected set of orbital 
elements were derived when considering the perturbational effects of earth oblateness 
and tangential atmospheric drag. This section analytically solves these equations by 
the method of two-variable asymptotic series. (To date, complete solutions have been 
obtained for both oblateness (J^ and J^) and oblateness/drag combined. ) Included are 
synopses of the FORMAC computer program used in obtaining the analytical solutions 
and the FORTRAN program used in numerically evaluating these solutions. 

4.1 GENERAL 

In the method classically known as general perturbations, six first-order equa- 
tions of motion can be formulated as functions of some fundamental jet of orbital 
elements. The perturbation effects are expressed analytically, and the element solu- 
tions are generally obtained by analytical integration of series expansions in one form 
or another. These solutions are explicit functions of time, constants of the problem 
and constants of integration. They define the vehicle state at any instant in time, as 
the epoch state conditions make the problem completely determinant. 

The primary difficulty in the general perturbation method has always been the 
overwhelming amount of analytical labor required to obtain the solutions. However, 
the state of the art in computer technology is such that automated manipulation languages, 
i. e. , languages for doing symbolic as opposed to strictly numerical mathematics, are 
now generally available. Consequently, many of these burdensome analytical tasks, 
such as series manipulations, function expansion, differentiation and integration, can 
now be alleviated. 

The language selected for use in this development is FORMAC (FOR MULA 
MA NIPULATION COMPILER). This language, currently available through the MSFC 
Computation Laboratory, was developed by IBM, and contains a wide range of ana- 
lytical capabilities (Reference 11). Consequently, it has proven itself a valuable tool 
for the application at hand. 
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4. 2 BASIC THEORY OF THE TWO- VARIABLE ASYMPTOTIC SERIES EXPANSION 
METHOD 


As was indicated earlier, general perturbation techniques employ series expan- 
sions for assumed element solutions. These expansions result in correspondingly 
expanded differential equations which are then analytically integrated. The type expan- 
sions employed in this study are classically known as asymptotic series expansions. 

It is the purpose of this section to provide an outline of the theoretical basis for such 
expansions, illustrating those concepts required in the particular application at hand. 

The discussion begins with some basic definitions and nomenclature. 


Definition 1 

Let f(t, €) and g(c) be real-valued functions, where £ is a small positive param- 
eter andj: ranges continuously over some set S of nonnegative reals. Then, a measure 
of the relative magnitudes of f(t, C) and g(c) may be obtained if a real (finite) K exists 
such that: 


Li m 
€“■> 0 


g(f) 


s K 


for all£ in S. Symbolically, the existence of this limit is denoted by writing: 

f(t, C) = 0(g(C)) 

which reads "f(t, €) is of the order of g(C). " The existence of the limit for all £in S 
makes this relation uniform in that K can be chosen independently of £. The function 
g(€) is called the gauge function, and when K = 1, f(t, C) is said to be asymptotically 
equal to g(€). If £ is a function of several real variables, the relation is said to be 
multivariable (Reference 12, pp. 180-185; Reference 13, pp. 1-3; and Reference 1, 
pp. 16-17). 

For purposes of clarification, consider the following example: 

Lcl t >o, o < c « 1 , f(t, e) = e 2 sin t and g(e) = e 
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Then clearly, 


f(t, e ) = 0(g(c)) 


since 


Lim 
f — *-0 


sin t 


* e 


for all t > 0; i. e. , uniformly in t_ as c— ►O . 


Another symbol, which is often used to measure the relative magnitude of two 
functions, bears a simple relationship to the order symbol of Definition 1. If f(t, e), 
g(e) and £ are as previously defined, then this alternate measurement is obtained 
when: 


Lim 

G 


iiL f) 

g(0 


= o 


for all t_ in S . Symbolically: 

f(t, € ) = o(g(e)) 

and is read "f(t, e) is small o of g(c). ” (When both symbols are employed, f(t, e) = 
0(g(e )) is often read "f(t, e ) is large O of g(e). ") 

The symbol small o, though not employed herein, is related to the large O of 
Definition 1 by: 


o(0(g(e))) = o(g(e)) 


Definition 2 


Let g^e ), i = 0, 1, 2, . . . , be a sequence of real-valued functions of the small 
(positive) parameter Then, this sequence is called an asymptotic sequence for 
e-^0 if, for each i_ (Reference 12, pp. 182-183; Reference 13, pp. 2-3; and Refer- 
ence 1, p. 17): 


Lim 


g.(0 


= o 


i 

Such a sequence is illustrated in the following example: 
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be defined by 


Let the sequence g.(e), i - 0, 1, 2, . . . , 


g ; ( 0 = c , X j+ i>\. >° 


for ail i_. Then this sequence is an asymptotic sequence, since 

g i4 1 (t) C Xi+1 

Lim r— = Li m — ; = 0 


-I * *41 . *-* * *** \ 

g.(0 A. 

-^o 1 e— 0 p 1 


for each i . 


Definition 3 

Let g.(e) and f^ (t) be real-valued functions of the small parameter e and the 
real nonnegative variable^, respectively. Then, the sequence of partial sums: 

N 


£ g <e)f f,) (t) 
i=0 


is called an asymptotic expansion to N terms of a function x(t, e) as e— «►() when: 

N 

X(t, e) = £ g.(e)i ' (t) + 0 (r. +1 (0) 


as c -►0. The asymptotic expansion is said to be uniformly valid when it holds for 
all t in some set S of nonnegative reals, i. e. , when 0(g (e)) is uniform in t . If 

_ _ i+i — 

t_ is expressed, at least formally, as a function of several variables, then the ex- 
pansion is said to be a multivariable asymptotic expansion. Such an expansion would 
have the form: 

N m-~ 

X(t,e) = £ g.(€)f w (t,t, ...) + 0(g. +i (f)\ 

as c-^0. For purposes of preserving the uniform validity of the expansion (Reference 

— rsj 

13, pp. 79-82, and Reference 1, p. 17), the variables t, t, . . . are taken as functions 
of c multiplied linearly by t. Here, F is termed the fast variable while t_ is termed 
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the slow variable (Reference 1, pp. 16-17). The first, second, and third approxima- 
tions to x(t, € ) are given as: 


g (C)f ,0) (t) + g (f)f (1) (t) 

0 1 

g 0 <f)f <0) (t> + gjff )f (1) (t) + g 2 (Of (2) (t) 

€ t 

To fix these ideas, take t > 0 and x(t, c) = e , where je is a small positive 
parameter. Let g^(c) = eVil and f^(t) = t 1 , i = 0, 1, 2, . . . . Then 


Lim 
C — *0 



Lim 
e— 0 



so that g.(€), i = 0, 1, 2 , ... is an asymptotic sequence. The sequence of partial 
sums: 

N 

is an asymptotic expansion to N terms of x(t, e) = e f *. 

Note that in this example the asymptotic expansion was convergent. However, 
there is to be no convergence requirement imposed on such expansions, and some ex- 
pansions may converge for some range of e, or may diverge for all £ . The practical 
applicability of the method is not determined by convergence of the series when 
i— •"», but by its asymptotic properties for a fixed value if i_ when e — ►O (Reference 
14, pp. 40-41). 



Hence, an important characteristic of asymptotic expansions is that the error 
made in approximating the given function by such an expansion is of the order of the 
first neglected term (Reference 1, p. 17). For this reason, it is important that one 
make a wise choice for the small parameter e when using this method. 
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Consider x(t, c) a function that is to be approximated by a two-variable asymptotic 
series expansion. Then, £ will be functionally related to two variables, say £ and £, 
in a linear fashion through _g_. Here again, £ will be termed a fast variable and £ a 
slow variable. Further, suppose that x(t, e) represents the solution to a differential 
equation whose independent variable is t_. To apply the technique, the differential equa- 
tion must be expressed as a function of both t_ (at least implicitly) and 

Thus, the initial value problem for an ordinary differential equation is converted, 
through use of a two- variable expansion, to one involving partial differential equations 
in £ and £. The two-variable asymptotic solution of the transformed problem will 
then involve certain undetermined functions which are defined by postulating that the 
problem possess a consistent asymptotic solution which is uniformly valid (at least to 
value of t_ of the order of the reciprocal of the small parameter). 

There are two concepts that aid in arriving at uniformly valid solutions, as 
opposed to those which are initially valid (i. e. , valid over some initial portions of 
their ranges). These are called the first and second uniformity conditions, respec- 
tively (Reference 1, p. 18). 

The first uniformity condition states that a multivariable asymptotic solution to 
a small parameter dependent differential equation cannot contain secular terms in the 
fast variable £ (i. e. , terms proportional to £), if the solution when e = 0 does not 
contain such terms. In short, if the solution to the differential equation when e = 0 is 
bounded in the fast variable, the solution procedure cannot unbound the solution when 
e / 0. Note that this condition is applicable only if the c = 0 solution is initially 
bounded (Reference 1, p. 18). 

The second uniformity condition is a result of the uniform validity requirement, 
and this condition states that: 

g 1+1 (c)f (l+1) (t, t) 

Lim — — = 0 

«— 0 g. (<)i (t, t) 

for each i_ and all t_ of some set S of nonnegative reals. Simply stated, the ratio: 
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cannot contain terms secular in the slow variable t_. This condition may be employed 
to eliminate nonuniform results even when the first uniformity condition cannot be 
applied (Reference 1, p. 18 and Reference 15, pp. 206-224). 

In Paragraph 4. 3, the two-variable asymptotic series expansion method will be 
employed in obtaining solutions to the variation-of-parameters equations derived in 
Section 2. Thus, the function x(t, 0 to he approximated by these two-variable expan- 
sions will represent some osculating element; § will be a small parameter arising 
through the pertuibational effects, and t, Twill be two time-scale variables associated 
with the time t . 


4. 3 APPLICATION OF THE TWO- VARIABLE ASYMPTOTIC SERIES METHOD TO 
THE DIFFERENTIAL EQUATIONS OF MOTION 

This method assumes that the solutions to the equations of motion can be ex- 
pressed as asymptotic series in two variables (t and t), i. e. , 


B = B (0) (t, t) 


+ e 


( 1 ) - ~ 
B w (t, t) + 


2 ( 2 ) - ~ 
e B' ; (t, t) + 


/ 0 ) - ( 1 ) - ~ 2 ( 2 ) - ~ 
e = e ' (t, t, ' + e e ' ' (t,t) + ( e ' (t, t) + . . . , etc. 

where B^, B^\ B^ 2 \ . . . , e^, e^\ . . . are functions of time (i. e. , solutions) as 
yet to be determined and: 


t = t (l+of 2 c 2 ) 

(fast variable) 

(4-1) 

T= €t 

(slow variable) 

(4-2) 


with o fg being an undetermined constant. 

In the asymptotic series expansion for a given element, the first term is re- 
ferred to herein as the first approximation to the total solution, and the sum of the 
first and second terms is referred to as the second approximation. For the element 


e, as an example: 


e = e<°> 

(first approximation) 

e - e<°> ♦ c e' l> 

(second approximation) 


1 — super-one solution 
— suner-zero solution 


These approximations will now be derived for the set of elements (B, e, i, Q, 
oj and M). First approximations will be obtained when considering both oblateness and 
oblateness/drag. Second approximations will be obtained in terms of the super-one 
solutions due to oblateness only, as it will be shown that the super-one solutions due to 
drag are negligible. The general procedure for obtaining the third approximations will 
be outlined. 

4. 3. 1 Obtaining the First Approximations to the Solutions 

The desired first (and second) approximations are obtained by solving the 

variation-of-parameiers equations (oblateness only or oblateness/drag) when consider- 

2 

ing only terms of the order of £ (i. e. , neglecting terms ). Since these equations 
are highly coupled, their solutions must be obtained simultaneously (at least in theory). 
However, by making reasonable assumptions, the solutions for each element can be 
obtained separately up to a point - this point being the formulation of a set of first- 
order ordinary differential equations having £ as the independent variable. To illus- 
trate the procedure leading to this point, the equation for a representative element 
will be considered in detail. 

4 . 3. 1. 1 Oblateness Only 

The element i_ is taken as the representative element, so it is necessary to 
expand each element appearing in Equation (2-36) to the first-order of £ . From 
Appendix B: 


4-8 


( 0 ) 

cos i = cost + € 


-i (1 > sini<°> 


+ € 


2 2 ( 0 ) 
sin u = sin u ; + c 


«< 1> sln2u (0 >| + < 2 [ ]♦... 


(l + e cos v)= (l + e^ cos i^) + € ] + £ 2 ]+..., etc. 

Therefore, to the first order of £, Equation (2-36) becomes, when considering 

only J_ and , (the solution procedure has not yet been extended to higher harmonics): 
2 J 


& = - € 
dt 


yiTr e 2 B (0) (l 


7 ' + e<°> cos 


W\ 7/2 


+ c 


2 (l-e ( °H 

r- 3 (0) 9 / J3 \ / (0) (0)\‘ 

# r e B ( 3?)( 1 + e 008 v ) 


S in21 (0) sin2u (0> 


(4-3) 




9/2 


.(0) (0) . 2 .(0) , 2 (0)\ 
cos i v cos u' (1-5 sin v sin u v '/ 


The solution method begins by assuming that Equation (4-3) has the asymptotic 
series solution: 



- ~ - ~ . (1) 

i(t, t) = i w (t, t)+*i (t, t) + ... 

(4-4) 

where 

? =t(i+*' a «*) 

(4-5) 


t = ct 

(4-6) 


Differentiating Equation (4-4) with respect to time yields: 


di _ gi dt | Q\ dF 
dt “ dt dt dT dt 


,( 0 ) 


:(1> 


/d£Z . . /< 

= (*t Sr~ •••/ dt v 


0i 

0t 


i(0) 


+ e 


<9i_ 

0? 


• ( 1 ) 


which becomes upon differentiating Equations (4-5) and (4-6): 


di /*i (0) .*£> \ , . 2) // W (1) 

df*Ur + , »r + -)( 1+ V 
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Rearranging in ascending powers of j[ yields: 


di 

dt 


3i^ + 

at + 


at a? 


1 2 

' 

+ 

-■*1 

. 


(4-7) 


Equating coefficients of like powers of c from Equations (4-3) and (4-7) results 
in the following partial differential equations: 



The problem has now reduced to solving these partial differential equations. 
Equation (4-8) implies that is either constant or a function of £ only. Conse- 
quently, 



di 


AO) 


dt 


= function of t or constant 


(4-10) 


In light of Equations (4-8) and (4-10), Equation (4-9) can be reduced to an ordi- 
nary differential equation if the constant of integration is considered a function of T. 
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The resulting solution in "integral form" is: 


,< 1 > 


d^ - frb r c 2 B (0>7 q ♦ e (0> c. - An 2 i (0) sin 2 u ,0) dt (4-11, 

■'*” 1 2 , 1 - e'< 2 




9 / T3\ /♦ 

V°> fe)« + e (0) cos ,< 0) >cos i (0) cos u ,0, a-5 sin 2 i ( °’sln 2 u (0) , dt 

' (0) 2a/z 
a - e l ' ) 


+ c (t) 

where C(t) is the integration constant. Before proceeding to solve the above integrals, 
it is desirable to transform the variable of integration from t to v This trans- 
formation is taken to be the standard Keplerian transformation (Reference 6, p. 221). 


2 3 2 

dy _ Jjl(l + e cos u) _ sTu B (1 + e cosM 

dt " 3/2 ” , 2 3/2 

P (1 - e ) 


Thus, 


-r— ^ 


2 3/2 


dv 


( 0 ) 


♦ e <0 »cos J°f 


(4-12) 


Substituting Equation (4-12) into Equation (4-11) und simplifying yields 

4 

(1) di ^ - r e 2 / B*°* (0) (0) (0) (0) , (0) 

i' = - — — t - I — “ (1 + e* ; cos 1/ ) sin 2 i' 1 sin 2u^ dis 


dt 


( 0) 2 2 
(1 -e' ’ ) 


+ r e 


©/ 


B 


(Of 


a - e (0 » , 


—3- (1 + e^cos cos i (0) cos 


(4-13) 


2 ( 0 ) 2 ( 0 ) ( 0 ) ~ 
(l-5sin i' WV ') di/ ' + C (t) 
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I 


Iri order to perform the indicated integration, it is necessary to know the 
dependence of the element functions B^, e^, i^ and (recall that = 
o,< 0) + v^*) upon However, these element functions are not yet known - in 

fact, the determination of these functions is the goal of the present development. 
Therefore, in order to proceed with the solution development, it is necessary to 
make a simplifying assumption based on the knowledge that the elements B, e_, i_ 
and vary slowly with time as compared to the element u. Specifically, it will 
be assumed that with respect to a dy^ integration, the element functions B^, 
e (0) . i^ and are constant. The effect of this assumption on the accuracj’ 
of the resultant solutions can be minimized by periodically rectifying the orbit 
and updating the epoch values of the elements, (As discussed in Paragraph 4. 5, 
an "updating procedure" is used when numerically evaluating the solution 
equations. ) 

hi "partial consideration" of this assumption. Equation (4-13) can be 
written as: 


j-CO „(0) 4 2 
.(1 di' * — B v r e 
i = t + 


dt 


( 0) 2 2 
2(1 V ' ) 


/ 


- (1 + e (0) cos 2 i (0) sin 2 u<°> ^ (4-14, 



A 


2 cos i^ cos u*°\l-5 sinV°*sin 2 u*%dl/ 0 ^ 


+ C (t) 


or in notational form as: 



“ t + K 2 (i) I 2 (i) + K 3 (i) I 3 (i) + C (t) 


(4-15) 


1 

1 

1 
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where 


V> * 


B ( 0 ) 4 r . 2 _ 

2(1 - e (0)2 ) 


(4-16) 


and 



(4-17) 



(1 + e (0> oos J°>, sin 2 i<°> s i» 2 u ( °W°> 


(4-18) 




T /-v , ,, (0) (0) 2 .(0) (0),, _ . 2.(0) . 2 (OK . (0) 

I 0 ( 1 ) = / (1 + e cos v ) cos i cos u (1 - 5 sin i sin u ) dp 


(4-19) 


(NOTE: In the above notation, the subscript on K_ and I_ indicates the earth-harmonic 
under consideration; the parenthetical (i) indicates the element i^ ) 


Since e^, i^, and are considered constant by the previously- stated 
assumption, inspection of the integrals given by Equations (4-18) and (4-19) reveals 
tnat each integral can be expanded to a series of single-term integrals of the general 


form: 


f(e 


(0) .(0) (OK 
1 , cc ) 


/■ 


. P 
sin v 


(0) Q tJ (°). (0) 

cos if dp 


(P, Q = 0, 1, 2, . . . ) 


which is directly integrable by "textbook” formulas. Unfortunately, such an 
expansion procedure results in many single-term integrals; to solve these by 
hand for each of the six elements would be an overwhelmingly laborious task. 
However, by utilizing the p.itomated tecnn jues of the FORMAC language, a computer 
pregram was written for the IBM 7094 to expand expressions similar to Equations 
(4-18) and (4-19) and then "solve" the single-term integrals by an identification and 


4-13 



substitution procedure. Basically, the program identifies through an iteration 
process, the values c i the exponents Pand Q occurring in each single-term integrand; 
it then substitutes the precoded solution for that particular integral. 

In general, the integrated solutions to Equations (4-18) and (4-19) will consist 
of terms secular in the independent variable and terms non-secular in 
i. e. , 

I 2 (i)=S 2 (i > I/ (0) + N 2 (l) (4-20) 

I 3 (i) = S 3 (i) l/ (0) + N 3 ( i) (4-21) 

where S denotes the secular terms and N the non-secular terms. The FORMAC 
program prints the answer arrays I_ , I , and S , S for each element; since 

m O Lt u 

these arrays are very lengthy, they are presented in Appendix C. 

In view of Equations (4-20) and (4-21), Equation (4-15) becomes 

‘ (1) = * ‘ + K 2 (i> [ S 2 (i)l/<0) + N 2 (i) ] + K 3 (i) [ S 3 (i)y<0) + N 3 (i) ] + C (?) 

(4-22) 

As shown in Appendix D, the element function v ^ is secularly related to the 
fast time-variable t^by 


( 0 ) ( 0 )_ 
v = n' t 


(4-23) 


Hence, the resolution of into secular and non-secular parts yields 


0 (0 0 0 - 0) 
v ' ' = v ' + t/' = n t + v ' 

s N N 


(4-24) 
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where is the non-secular part of vet to be determined. 

N 

Substituting Equation (4-24) into (4-22) yields 


i' 1 * = - — t + K (i) S (i)n ,0) t + K (i) 
dt * 


i)[ s 2 ( ,) . N <0) + N 2 ( i) ] 


+ K 3 (l) S 3 (i)n (0, t + K 3 (i) 


S 3 W 


l > V N <0) + N 3 (i) . 


+ C (t) 


which becomes after rearranging 


i (1) = 


di 


: ( 0 ) 


dt 


+ K 2 (i) s 2 (i) n (0) + K 3 (l) S 3 (i) n <0) 


(4-25) 


+ K 2 (i) [ S 2 (i)l 'N < ° ) + N 2 (,) ] + K 3 (i) [ S 3 (i> ^ + N 3 (i) ] + C <*> 

At this point, the first uniformity condition (see Paragraph 4.2) can be imposed. 
Essentially, this condition requires that any approximate solution to the element 
i_ not contain a secular term in the fast variable T since the solution to the 
differential equation for i_ did not contain a secular term when 0. In order 
for this condition to be satisfied, it must be that 


^T~ + K (i) S (i)n (0) + K (i) S (i)n (0) 
dt 22 d 3 


(4-26) 


In view of Equation (4-26), Equation (4-25) becomes merely 


i (1) = K 2 (i) 




(i) 

o N 


( 0 ) 


N 2 (i) 


«] + K 3 (1) [ 




( 0 ) 


N 3 (i) 


*>] 


+ c (t) 


(4-27) 
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Now, as is evident in Appendix D, one method for obtaining the non-secular 
part of (i.e., v ^ ) would be to evaluate the indicated Fourier series 
This is not necessary, however, since Equation (4-24) can be rearranged as 


( 0 ) ( 0 ) 


N 


= V 


n (0) T 
- n t 


and and n (0) T are known. 


Thus, Equation (4-27) can be written as 


i* 1 ’ = K 2 (i) 


1)^(1) v (0) + N 2 (i)J - K 2 (i) S 2 (i)n <0) t + K 3 (i)rS 3 (i)i/ (0) + N g (i)l 


( 0 )- ~ 
-Kg (I) S 3 (i)n l T+C (t) 


which becomes by Equations (4-20) and (4-21) 


i (1) - k 2 (1) 


I 2 (i> - S 2 (i)n ,0 ¥ 


+ K 3 (i) 


i) I 


3 (i) - S 3 (i)n ,0) tl + C (t) 


(4-28) 


where K (i) and K (i) are given by Equations (4-16) and (4-17), and I (i), I (i), S (i), 

and S g (i) are obtained from the FORMAC program (see Appendix C). It should 

be noted that although the appearance oft_ in Equation (4-28) suggests secularity, 

this secularity is "cancelled" by that appearing in I (i) and I (i). Consequently, 
rt) - 2 3“ 

is non secular in t, thereby satisfying the first uniformity condition. 

Returning to Equation (4-26), it follows that 


= K 2 (i) S 2 (i)n (0) + K 3 (i) S 3 (!)n <0> 
dt ° 


From the FORMAC results presented in Appendix C 


s 2 (i) = 0 

S 3 (i) = e (0) cosi (0) cos (1 -|sin 2 i<°>) 


(4-29) 


(4-30) 

(4-31) 


4-16 


Substituting Equations (4-16), (4-17), (4-30) and (4-31) into Equation (4-29) yields 


di_ 

dt 


:( 0 ) 


B 




3 ( 0 ) ( 0 ) .( 0 ) ( 0 ) 

ife n e cos i' cos a; 5.2 (0) 

73 a “ 7 i ) 

(1 - e (0)2 , 


This is the first-order ordinary differential equation having T as the independent 

variable which was referred to at the beginning of this section. Using a procedure 

identical to that just illustrated for the element i, the corresponding equations 

for the remaining elements can be formulated. The set of equations for all elements 

is presented below, along with the (approximate) solutions to the equations as 

derived in Appendix E. These solutions were obtained by a method set forth in 

Reference 16, whereby e and oj_ are considered to vary simultaneously and terms 
2 

of the order of e_ (or smaller) are ignored. The constants (A, a, C , C , . . . , C ) 
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appearing in the solutions are defined in Appendix E. 


Element e 


de<°> 

dt 


B<°> 6 r 3 n (0) cosw (0) sini (0) 

W (5 s in 2 i«»-l 

( 1 _ e (0)2 ) 2 \4 I 


(4-32) 


(' i c 1/2 

e (0) = ^A 2 + 2 A sin (C 2 t+ a) +| _L| J 


(4-33) 


Element co 


d J0) B(0) 4 rp 2 n (Q) 6 

dt = (l- e (0)2)2 (2— 2sin 2 i (0) ) + 


" (0)6 (if) -e 3 n"" 


(1 - e (°)2)3 

- e(°) esc i(°> sin + -L- sin i«» sin J°) (1 -5. sin 2 i«»)l 
e (0) 4 J 


[%- e(°) sini(°)sin w«» cos 2 i<°> 


«(°)* tan-1 


c 

A sin 'C 2 t+a) + —L 


A cos (C 2 t+a) 


(4-34) 


(4-35) 


Element i 


. (0) B<0) 6 |^a _) r e 3n(0)e(°)cosi(0) cosaJ (0) 


d 

dt 


( 1 — e * 0 ) 2 ) 3 


( 1-1 sin 2 j'0)| (4-36) 
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( 4 - 37 ) 


( 0 ) 

i 


.(0) 

1 0 ' 


' ( 0 ) . 
e sin u> 


( 0 ) 


( 0 ) 

— e q sin u> 


(0)\ 

°) 


Element 0 


(0) 4 2 (0) (0) (0) 6 /43 


d n 


(0) B r e £ n cos i 


B <0 ’ r e ^ n<°>e«»sin W «» 


d7 (1— e (0)2 ) 2 




(ctn i (0 > - sin 2 i <0) ) 


o (0) _<0) C 5 / ,0) (0) (0) (0>\ / q \/~ -v 

“ = n 0~C^\ e cosw ~ e 0 cosw 0 ) (°5 C4 )( o) 


Element B 


dB<° > 

d"t 


= 0 


NOTE: Since 


B<°> = B ( J> 


n = y/n 


it follows tnat 


( 0 ) 

n 




(Or 


— Tjf 0 ' 

- VM B Q 


(0) 
= n 0 


Element M 


dM<°) 
d t 


B (0) 4 r 2 (0) 4B(0)6 (7 i )r e 3 n(°> 

= 15— T /j_l .2.(0)) \ J 2/ 6 

U-e(0)2)3/2 U 2 ' 3e<0) (l -^2) 5/2 

+ | sin i (0 > gin to <0 > ( 1- | sin 2 i<°»)J 


j-3e(°' 2 sin i* 0) sin (1--| sin 2 i<°) 


( 4 - 38 ) 

( 4 - 39 ) 

( 4 - 40 ) 

( 4 - 41 ) 

( 4 - 42 ) 

) 

( 4 - 43 ) 


4-18 


(4-44) 


M"'’ = [n%- V ] • M,/ 0 ' '-P- -4^-MCt-c. 


C’,C c\c„ 

1 y 1 ” 


2 


(C K ' 4< V (0) (0) <(>) (0), 

; (c cos u." - c cos a.' ) 


NOTE: The first term on the right-hand side of Equation (4-44) is the Keplerian 
change in M_ that takes place during the time interval (t—t^). When applying the 
asymptotic series solution method to Equation (2-39) to obtain Equation (4-43), 
the Keplerian variation is ignored since this variation can be solved in a straight- 
forward manner from the Keplerian equation 


dM 

dt 


= n 


Consequently, the Keplerian change must be added to the solution of Equation (4-43). 


Element v 


As mentioned in Paragraph 2.1, the element v is obtained by a Fourier-Bessel 

2 

expansion involving M and e. To the order of e_, this expansion is (Reference 6, 
p. 89) 


(0) (01 (0) (0) 5 (0) (0) 

t/ = M' ' + 2e* ; sin ’ + -e' ’ sin 2M l ; (4-45) 


4. 3. 1.2 Oblateness and Drag 

The element e is taken as the representative element, so it is first necessary 
to form the composite differential equation for e when considering oblateness and 
drag. Since perturbative forces are additive, this is done by merely adding 
Equation (2-35), J and J terms only, to Equation (2-42). It is then necessary to 

M o 

expand each element in the composite equation to the first order of e. Using the 
expansions presented in Appendix B results in 
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i 


de 

dt 


V ^ rB (°) 7 r e 2 (l + f «> lc O , V ( 0 >) 3 r 

’ e (l _ e (0)2)7/2 [ iin ‘' < ° ) (1 ♦ *<°> cos ^<°>) ( 1 - 3 sm 2 i<°> ,m 2 u<°>) 

♦ tin 2 sin 2 u*°* ( 2 coi i4°) + e*®l cos 2 y<°l +e^*)J 

4 > Bl 0 ) 9 (^) r e 3 ( 1 + e ( 0 > co,^»f 

~ * ' (l-«.«>>2y9/2 


Jsin i ,0) sin u ,0) sin v <0) (l+« ,0) cotv (0 *)f3 — 6 sin 2 i (0> *in 2 u ,0) ) 


- f ( 1 - 6 sin 2 i<°> sin 2 u<°>) ( 2 cos k‘°> + e«» cos 2 ♦ e <°>) (cos u<°> sin 

K* p B<°> 4 (e ( °) ♦ cos + e (0) 2 + 2e<°> cos p<°>) 1/2 

* ‘ n (0> (1 _ e (0)2 ) l/2 


The solution method begins by assuming that Equation (4-46) has the asymptotic 
series solution 

e (t, t) * (t, t) + ee^ (t, t) + . . . 


Following the procedure outlined in the previous section for the element i the 
solution in ’’integral form” is obtained (corresponding to Equation (4-14): 


eUl.-de^’ j 


B ( 0) 4 .2 


dt 


e (0|2j2 


/ (1+ c<°> cos v<°>) |sin 1 * e (°> cos v <°>)(1 - 3 sin 2 i'°> sin 2 u<°>) 


+ sir 2 i<°) sin 2 u<°> (2 cos V (°) + e«» cos 2 ^°> ♦ e(°>)J d^°> 

3( I — f'°*2 j3 " / f 1 ’ cot jam | ,0 > gin u (0) tin e 101 (1 + l - osp ,0> ) 

( 3 — S sin 2 i*^ sin 2 


- | (l -5sin 2 i< 0 >,m 2 u<°>)(2coss'< 0 » + e (0) co ,2 v (0) < . e (0,| C os u<°> s.n i<°>] dK<°> 
+ K»(l— e<°> 2 ) P(e<°>^cosi/< 0 >)(1 ♦ e<°> 2 » 2 ><°> cos ,«») 1/2 

^B<°>2 J (1 ♦e<°><W°» 2 dt’ (0) + c('t) 


(4-46) 


(4-47) 

I 

1 
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or in notational form as 

/i) 

* =_ dT~ 1 +K 2 (e) ! 2 (e) + K 3 (e) ¥ e > + K D (e) l D& + C (t) (4-48) 

where 


B «» 4 ,2 

K 2 (e) (l— e(0)2}2 

(4-49) 

4B(0,6 (j) r ‘ 3 


K3(e,= auii 

(4-50) 

Y (r , . K*(l-e<°) 2 ) 

(4-51) 

l2( e ) = / — (l + cos i>^)|sin J^°l(l + cos (1 — 3 sin 2 j(°) sin 2 u^ 0 ^) 

(4-52 

+ sin 2 i*°) sin 2u^) (2 cos y(°) + cos 2 + e^))J d^ 0 ) 



13 (e) (1 + e^°) cos 2 |sin i ^ sin sin + e* 0 * cos ^°^)(3 — 5 sin 2 sin 2 i/®)) 

~ ^ (l — 5 sin 2 i^) sin 2 u^) (2 cos cos 2 + e^) cos sin i^lj (4-53) 


f p (e^) + cos (l + e(®) 2 + 2e^) cos I'*®)) ^ 2 <n , 

W” = 1 »^™7,<0>,2 ' 


(4-54) 


(NOTE: In the above notation, the numerical subscript on j< and I_ indicates the 
earth -harmonic under consideration, the subscript D indicates drag, and the 
parenthetical (e) indicates the element e. ) 
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Since e^, i and (recall that u^ = o/°^ + t/% are considered constant 
by the previously stated assumption, inspection of the integrals given by Equations 
(4-52) and (4-53) reveals that each integral can be expanded to a series of single- 
term integrals of the general form 


,, (0) .(0) (ok r P 
f(e , i v u? ) I sin v 


( 0 ) 


Q„(°h (0) 

cos v d 


(P. Q = 0, 1, 2, 


which is directly integrable by "textbook" formulas and, therefore, the FORMAC 
program. However, in order to readily integrate Equation (4-54), it is convenient 
to employ the binomial series approximations 


,, ( 0 ) ( 0 ) ( 0 ) 1/2 

(1 + e + 2e cos v ) 


„ ( 0 ) ( 0 ) -2 
(1 + e cos v ) 


( 0 ) ( 0 ) , ( 0 ) 2 ( 0 ) 
1 + e w cos v 7 + 1/2 e' sin u 


( 0 ) ( 0 ) ( 0 ) 2 ( 0 ) 

1 - 2e cosy + 3e cos v + 


+. . . (4-55) 
. . . . (4-56) 


It is also necessary to know the functional dependency of atmospheric density_£ 
upon true anomaly In the past, this dependency has been established by 
using very «imple models of atmospheric density, such as an exponential model or 
a power-law model. Though convenient to work with, these types of models do not 
provide realistic simulations of the actual environment since they are structured 
to represent the variation of density with altitude only. Density actually varies 
with solar and geomagnetic activity, time of year and position relative to the 
sub-solar point (diurnal bulge), as well as with altitude. 

Realistic simulations of long-term satellite motion must include these additional 
variations in the density model. For example, using a simple density model (the 
1959 ARDC) to compute the lifetime of Satellite 1961c results in a lifetime of 179. 1 
days. The actual lifetime was 525.5 days - an error of 66%I On the other hand, 
using a realistic model (the 1970 Jacchia) produced a lifetime of 537.9 days; an 
error of only 2.4%. 

The difficulty with using a realistic density model is in expressing density as 
a function of true anomaly. An examination of the 1970 Jacchia model shows 
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how complex a realistic model is and, consequently, how difficult it would be to 
implement directly into a general perturbation technique. Yet, in order for the 
general perturbation technique to be as accurate as numerical solutions, it is 
desirable to use the 1970 Jacchia model. 

A rather unique approach to the use of a realistic density model is taken in 
this study. Specifically, the variation of pwith is approximated by the Fourier 
series 

4 

p = 1/2 a Q cos ki/°* + b^ sin ki^ 0) ] (4-57) 

k=l 

where a^, a^, and are Fourier coefficients determined in the following manner: 

A table of density values is computed for intervals of true anomaly around one 
orbital revolution by numerically evaluating the 1970 Jacchia model. Integrals 
associated with determination of the Fourier coefficients are then computed by the 
Trapezoidal Rule. (It was found that the Fourier series using coefficients 
through^ and b^ give an excellent approximation to the functional dependency of 
density upon true anomaly. ) 

Because of the dynamic nature of the density function, the series approximation 
will not hold for long periods of time. (In fact, this is one area in which further 
study is recommended - see Section 6. ) The length of time depends somewhat upon 
the amount of resolution in the density input data (solar flux, geomagnetic index, 
etc. ) and upon the orbital conditions. For instance, if daily values of solar flux 
and heating parameters are used, the series would need to be evaluated at least 
daily. If the orbit is in a state of rapid decay, the series could require more 
frequent evaluation. As discussed in Paragraph 4.5, the Fourier coefficients are 
updated at required intervals when numerically evaluating the solution equations. 
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Returning to the solution procedure and substituting Equations (4-55) through 


(^-57) into Equation (4-54) yields 

I, 


w"/- 


(1/2 a + a cos v ^ + ...+ b, sin4 (e^ + cos t/% 

0 1 4 (4-58) 


(1 + e (0) cos „< 0) + l/2e (0) sinV 0) > (1 - 2e (0) cos „ (0) + 3e (0) cosV°>) 


dv 


( 0 ) 


J 

The FORMAC program is utilized to expand Equations (4-52), (4-53) and (4-58) 
and then "solve" the single-term integrals. In general, the integrated solutions 

will consist of terms secular in the independent variable and terms non-secular 

. ( 0 ) . 
in v ; i.e. , 


X 2 (e) = S 2 (e) ^ 0) + N 2 (e) 

(4-59) 

ye) - S 3 (e)y (0) + N 3 (e) 

(4-60) 

I (e) = S (e)l/°^ + N (e) 
D D D 

(4-61) 


where ^denotes the secular terms and N the non-secular terms. The FORMAC 
program prints the answer arrays L, L* In S , S for each element; 

since these arrays are very lengthy, they are presented in Appendix C. 

In view of Equations (4-59) through (4-61), Equation (4-48) becomes 




t + K 2 (e) 

S 2 (e) l / (0) + N 2 (e) 

+ K 3 (e) 

S 3 (e)l/ (0) + N s (e) 


. 


. , 


S ( e)u 
D 


( 0 ) 


+ N D (e) 


+ K D (e) 


J 


+ C (t) 


(4-62) 




As shown in Appendix D, the element function is secularly related to the 
fast time -variable t by: 

v 


(°)_ ( 0 ) ; 
„ = n' ; t 


(4-23) 


,( 0 ) ... 


Hence, the resolution of v' into secular and nonsecular parts yields: 

rot (°) (0) tot- (°) 

V W) = v s + v N = n(°) t + v N ' (4-24) 

where is the nonsecular part of yet to be determined. 

Substituting Equation (4-24) into Equation (4-62) and rearranging yields: 


,d) = 


--7T + K 2 ,e) S 2 (e) n<0) + K 3 (e > S 3 < e > n(0) + k D (*) Sd ,e > n(0) 


(4-63) 


dt 
+ K 2 (e) 

+ K d (e) 


(0) 

®2 ( e ) v N + ^2 + K 3 (e) 

_ (0) 

v N +N D( e ) \+C(l) 


S 3 (e) v 


( 0 ) 

N + N 3 (e) 


At this point, the first uniformity condition (see Paragraph 4.2) can be imposed. 
Essentially, this condition requires that any approximate solution to the element e 
not contain a secular term in the fast variable t, since the solution to the differential 
equation for e did not contain a secular term when e = 0. In order for this condition 
to be satisfied, it must be that: 

= 0 

(4-64) 


~*fr + K 2 (e) S 2 (e) n(0) + K 3 S 3 (e) n(0) + K D < e > S D ( e ) n(0) 


In view of Equation (4-64), Equation (4-63) merely becomes: 


e^> = K 2 (e) 

+ K d (e) 


(0) 


(0) 

s 2 (e) V N + N 2 (e) 

+ K 3 (e) 

S 3 (e) v N + No (e) 


(4-65) 


( 0 ) 

S D ( e ) v n + N d (e) 


+C(t) 
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Now, as is evident in Appendix D, one method for obtaining the nonsecular 
part of v^(i.e. , v^) would be to evaluate the indicated Fourier series. This is 
not necessary, however, since Equation (4-24) can be rearranged as: 


N 

and v<°> and n^ t are known. Thus, when considering Equations (4-59) through 
(4-61), Equation (4-65) can be written as: 


e (1) = K 2 (e) 

+ K D (e) 


I 2 ( e )-S 2 (e)n(°U 


+ K 3 (e) 


h (e) “ S 3 (e) n(°) t (4-66) 


n(0) t 


+ C(t) 


It should be noted that, although the appearance of t in Equation (4-66) suggests 

secularity, this secularity is cancelled by that appearing in I (e), I (e), and 

( 1 ) “ Z 3 
I D ( e ). Consequently, e is nonsecular in t_, thereby satisfying the first uniformity 

condition. 


Returning to Equation (4-64), it follows that: 

- K 2 (e) s 2 (e) n<°> + K 3 (e) S3 (e) n<°> + K D (e) S D (e) n<°) 

From the FORMAC results presented in Appendix C: 

S 2 (e) = 0 

S3 (e) = -|-sm cos JO) (1— e ^) 2 ) (Ij-sin^ i(®) — 1) 

S D (e) = +|b 3 ) + (~--a 0 +-La 2 ) e(°) 


(4-67) 


(4-68) 

(4-69) 

(4-70) 


(Recall that a , a , a, and b are the Fourier coefficients appearing in the density 

U 1 Z u 

function approximation. ) 
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Substituting Equations (4-49) through (4-51) and Equations (4-68) through 
(4-70) into Equation (4-67) yields: 


de^®) _ 


B<°> 


6 /±l\ 


r 3 n (0) s j n j(0) ^ JO) 


dt 


(1— e (0)2)2 
n<°) K* (l-e^°) 2 ) 




(|sin 2 i(0> - ij 

HS + f fc s) t ( _ T*0 + T* 2 j' 


,(«) 


This is the first-order ordinary differential equation having t_ as the independent 
variable which was referred to at the beg innin g of this section. Using a procedure 
identical to that just illustrated for the element e, the corresponding equations for 
the remaining elements can be formulated. The set of equations for all elements 
is presented below, along with the (approximate) solutions to the equations as derived 
in Appendix F. The constants (a, b, c, X X C , C . ... c_, D n , B , 3 , 6 , 

1 Z 1 Z 8201 —1 

and 5 ) appearing in the solutions are defined in Appendix F. 


Element e 


de<°) 

dt 


b(°) 


*&)-*■ 


sin i^) cos c*/®) 


where 


(l-e(°)2)2 

+ n (°)K*(l-e(°)2) 

H B' 0 ) 2 

e (0) = (f 2 + ,2)1/2 

. bt 

€ =6 xp [Xj cos c t + X2 sin c t] 
bt . 

7 ~ e X p [~ Xj sin c t + X2 cos c t ] — 


|-|sm 2 i<°> - lj 

( _ T a l + T b s) + (~T a ° + T 32 ) * <0) 


(4-71) 


(4-72) 


— ab 
b 2 + c 2 


ac 

2 a „2 


b * + c 
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Element w 


*gg> - B<0>4 r ^ 2 n<0) (2-5.taM»)) 

dt (i— e (0)2»2 \ 2 / 


B (0)6(i^r e 3 n (0) 


(l_ e (0)2 ) 2 \ 2 / (1— e <0>2>3 

— e^®) esc i^®) sin oj^®) + — — sin i^®) sin u^®) (l —Asm 2 i^®)) 
JO) \ 4 / 


Me(®)sini(®)sino;(®)cos 2 i(®) 

4 


+ n (®)K*q-e(®> 2 ) 

B(°) 2 e(®)/4 


fi (0) 

-l bl -3b 3+ lb 2 e(®) 
2 1 2 6 4 1 


(4-73) 


JO) = tan 




where .? and rj are given on the previous page. 

Element i 

g(0)® /£3\ r 3 n (0) e (0) cos i(0) cog w (0) 

dj(?L . \hl ll-|.in 2 i«») 


dt 


(l-e(®) 2 )3 


(4-74) 


(4-75) 


mi (®) c 3 i fOi ■ fO) (®) • (®>\ 

- i 0 +_i I e' u J sin CJ^ U, — e Q sin w 0 J 

Co ' 


Element q 


(4-76) 


,0,4 2,0, .,0) B(«) 6 (i),/n<») ( l0l- U (0) 

d j2(0) B<°> t e 2 n (0) cos i< u > + _ \*2l ___ 


dt 


( l_ e (®)2)2 


( l_e(®)2)3 


|ctn i^®) — sin 2 i^®) j 


(4-77) 


= — p 0) cos w ( ®) — e^coso/JJ* j + | c 5 C 4 j *0^ 


(4-78) 
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Element B 


dB<°) 

at 


■ n«»K> 
2B<°> m 



(4-79) 


B<°) 


( 0 ) 

* 0 D 2 


( 0 ) . . 

(@0 + e 0 )(t - t 0> + B 


<0> 


(4-80) 


NOTE: Since 


it follows that 


n “v^TB 3 


n(°)-v^ (0)3 


Element M 


(4-81) 


3e^ 0 ^ 2 sin sin u>^ ( 1 — — sin 2 i^ 0 *) 
4 


, 4B< 0 > 6 /-2\ r 3 n<°> r 

d M(0) ,B^^ (l _ ltin2i(0) \ 4B ( jjh jL f 

dt (1 _ e (0)2 } 3/2 \ 2 I 3e (°) (1 _ e (0)£)5/2 [ 

4* “ 01 * “' 0I (* -f “” 2 i<01 )] * ‘" K Sl g [ft 6 * *!Sto»*(^‘ *7 b2 ^ b3 ) 


(4-82) 


M<°>- [n<°> <t— to)] ♦M ( S ) + c 6 + (4C 7 -C 8 ) bT ^ + °2 6 0 + 6 1 *oj 

e bi I - - 1 

.. *J? — 1 (X 2 c— Xj b) sin c t ♦ (X 2 b + Xjc) cos c t , 

O 2 + c 2 I ’ 


[t-y 


— (4 C, — Cg) 


e bt 0 
_ e *P 

b 2 + c 2 


|(X 2 C — Xjb) sin c Iq + (X 2 b + X^c) cos c t j 


(4-83) 


NOTE: The first term on the right-hand side of Equation (4-83) is the Keplerian 
change in M that takes place during the time interval (t - t Q ). 
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Element v 


As mentioned in Paragraph 2.1 t the element V is obtained by a Fourier-Bessel 

2 

expansion involving M ana e. To the order of e_, expansion is (Reference 6, 
p. 39): 

„(0) = M (°) + 2e( 0 ) sin sin 2 M^) (4-84) 

4 


4.3.2. Obtaining the Second Approximations to the Solutions 

It so happens that the second approximations arc very nearly obtained during 
the process of deriving the first approximations, since the super-one solutions are 
merely functions of the super-zero solutions and an integration constant. The pro- 
cedure for completing the derivation of the second approximations will now be 
illustrated. As mentioned at the beginning of Paragraph 4.3, and as will be more 
thoroughly discussed at the end of this section, drag need not be considered since 
the super-one solutions due to drag are negligible. 


The element i_will again be considered in detail as a representative element. 
Recall Equation (4-28): 


i(l) =K 


2 (*) 


I 2 (i)-So(i)n<°>i 


I 

+ K 3 (i) 1 1 3 (i) — S 3 (i) n<°> t 


+ C(t) 


(4-28) 


where K (i) and K (i) are given by Equations (4-16) and (4-17), n^ is given by 

Z o 

Equation (4-81), and I (i), I (i) , and S (i) are obtained from the FORMAC program 

2 3 3 

(see Appendix C). Hence, once the first approximations are known, i can be 
computed from Equation (4-28) after the constant of integration C (t) has been deter- 
mined. In theory, a second application of the first uniformity condition (see Para- 

/v 

graph 4.2) would provide a means of determining C(t); unfortunately, this requires 
at least a partial formulation of the third approximation (see Paragraph 4.3.3). 
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To readily proceed with the solution development for i^, it is convenient to 
make a simplifying assumption based on the supposition that C(t) , being a function 
of the slow time-variable t_, is slowly varying itself. Specifically, it will be assumed 
that C(t) is a constant, C, The effect of this assumption can be minimized by periodi- 
cally rectifying the orbit and updating the epoch value of this constant. (As discussed 
in Paragraph 4. 5, an updating procedure is used when numerically evaluating the 
solution equations and, as shown in the plots of Appendix H, the C (t) for each element 
remains sufficiently constant over time intervals which are not extreme. Further- 
more, as discussed in Appendix G, these plots well describe the functional form of 
the C(t)'s obtainable when considering the third approximation. > 


In considering this assumption. Equation (4-28) can be written as: 


i (1) =K, (i) 


I 2 (i)-S 2 (i) n'°>t 


+ K 3 (i) 


I 3 (i)-S 3 (i)n<°>t 


+ C (4-85) 


The constant C can now be evaluated from initial (or epoch) conditions. From Equation 
(4-14), it can be seen that C is the constant associated with a dv^ integration in 
which all other element functions are considered constant. So, at epoch time t^. 


Equation (4-85) becomes: 

K 2 0) J h <*> ~ s 2 (*) n(0) 1 1 + K 3 <*> 1*3 (*> “ S 3 <»> n(0) 1 \ ! t A + C 


.(1) 

1 0 


I 




(ro 


(4-86) 


where [ ]_ 

*0 

( 0 ) 

indicates that the functions of v within the bracket are to be evaluated using the 
epoch value v^. Functions of the other elements (such as sin u/ 0 ^) are evaluated 
using current values. For example: 

[sin cos = sin cos 

e 0 
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(The FORTRAN program discussed in Paragraph 4. 5 employs a procedure for 
updating this epoch value. ) 


In view of Equation (4-86), Equation (4-85) becomes: 
i ( 1 ) -i ( 0 )+K 2 (i> I 2 (i)-S 2 (i)n( 0 >t + K 3 (i) I 3 (i) — S 3 (i) n(°> t 
K 2 (0 j I 2 0) ~ S 2 (i) n<°) t j + K 3 (i) J I 3 (i) - S 3 (i) n<°> t J 


(4-87) 


Since K and K 

& O 


are not functions of 



Equation (4-87) can be written as: 


iW-i^ + KjO) I 2 (i) - S 2 (i) n<°> I 


+ K 3 (i) I 3 (i) — S 3 (i) t 


- K 2 (i) |l 2 (i) - s 2 (i) n<°> tj - K 3 (i) |l 3 (i) - s 3 (i) „<0) { j _ 


or more concisely as: 


0 + k 2 (i) I 2 (i) ~ s 2 (i) n<°> t + K 3 (i) I 3 (i) - s 3 (i) n<°> 1 1 ! 

J t 0 


(4-88) 


where K g (i) and K 3 (i) are given by Equations (4-16) and (4-17), n* 0) is given by 
Equation (4-31), and I 2 (i), Ig(i), S 2 (i) and S 3 (i) are obtained from the FORMAC 
printout (see Appendix C). 

To this point, the solution for i^ has been considered in notational form. For 
a more revealing look into the actual solution, it is necessary to substitute Equations 
(4-16), (4-17), and (4-81) and the FORMAC results I 2 (i), I 3 (i), S 2 (i), and S 3 (i) into 
Equation (4-88). The solution resulting from these substitutions is presented in 
Equation (4-89). 

Using a procedure identical to that just illustrated for the element i_, the cor- 
responding equations for the remaining elements can be obtained. The equations for 
all six elements are summarized in rotational form following Equation (4-89). 
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Solution 


( 4 - 89 ) 


(1) (1) B<°> 4 r e 2 

l * J n ♦ £— 


2 ( 1 — *( 0 ) 2)2 


-2 e<°> sin 2j(°) sin v*°> sin w<°> cos w(°> + 4e ,0 > sin 2i ( °) sin 2 i>(°> sin 2 w<°> cos i><°> - 2 e (0 ) sin 2 i<°> sin 2 *<°> , 


„( 0 ) 


♦ le(°) sin 2 i(°> sin 3 u (0) ros JO) _ | e (0) sin 2 j<0) sin 2 JO) cos „(0) + | e (0) sin 2 j( 0) cos r <0) _ 2 sin 2 ,(0) ^ „<0) sin JO) ms JO) 

cos w<°> - sin 2 i(°> sin 2 *< 0 > ♦ 2 sin 2 i<°) sin 2 *«» sin 2 a>(°) ♦ 2e«» sin 2 i<°> sin ^ sin w <°> cos - i*<°> sin 2 i<°> sin 2 ^ sin 2 u <0) 

cos ^o 1 ♦ | e ( °) sin 2 i<°> sin 2 ^ cos ^ - i e <°> sin 2 i<0) ^3 „• J> sin w (0) ros w (0) + 4 e (0) sin 2 s <0) ^2 ^(0) cos W _ 2 e (0) ^ 2 *(0) 

J 3 ° 3 

cos u ♦ 2 sin 2 i<°> sin v JJ* sin w<°) cos /J* cos w<°> + sin 2 i< 01 sin 2 v JJ* - 2 sin 2 i<°> sin 2 j/JJ* sin 2 ui(°> | 

(j J ) r e 3 , 

2 -5 sin 2 i<°) sin p<°> sin 2 w< 0 > cos i«» cos w<°> - 5 sin 2 i<°> sin 2 ^°> sin cJ°> cos i<°> cos „«» + 20 sH 2 ^0) ^2 „(0) sijl 3 w (0) 


B (0) 6 / J 3\ , 3 
(1 —#( 0 ) 2)3 


cos i<°> cos p(0) + M sin 2 i< 0 ) sin 3 v™ sin 2 JO) cos i( 0> ^ w (0) _ | sin 2 i( 0) sin 3 „<0) cos 4 (0) ctr JO) _ 5 si „2 i( 0) sin 3 JO) m i<0 ) ^ JO) 

3 

♦ sin *<°> cos i(°) cos u/0) + sin o/°> cos i«» cos k«> ) -1 e (°) c«» sin 2 i(°» «>, i<0> ^(0) + e (0) „(0) ^ i( 0) ros _ 10 e (0) ^2 ,(0> sin ,,(0) 

sin 2 w(°) cos i(°) cos i-* 0 ) cos J°> + |c(°) sin 2 i<°) sin if*°> cos j(°) cos ^°) cos u>(°> - 10 c(°) sin 2 i<°) sin 2 *«» sin J°> cos i<°) 

+ 15 e(°> sin 2 i(°> an 2 v(°> sin 3 a,<°> cos i<°> + 10 e(°) sin 2 ««» sin 3 c<°> sin 2 c*/°> cos i(°) cos i/(°> cos u>(°> - £ *(0) sin 2 j(0) jin 3 „<0) 


cos i(°) cos V <°: cos W (°> + i£ e(°) sin 2 i(°> sin 4 *«» sin <J°) cos i<°) ♦ 10 e(°> sin 2 j(°) sin 4 v<°> sin 3 cj<°> cos i<°) * *«>) ^ j 


( 0 ) 


cos i(0) cos y(°) cos cj(°) — e(°) sin 2 i/(°) sin cj ( °) cos i( 0) 




Solution (continued) 


♦ 5 sin 2 «<°> un / sin 2 JO) ros *(0) m JO) + 5 ^2 j<0) ^2 „< jj> ^ JO) ^ j<0) ^ J jj> _ 2± ^2 *(0) ^2 ^3 JO) „„ 4 (0‘ ^ JO) 

-M sin* i«» ^3 ^2 JO) ro8 j( 0) ros JO) + | sin 2 j(0> sin 3 ^ ^ *(0) ^ W M>) ^2 j( 0) ^3 JO) ^ i( 0> w JO) _ sin JO) ^ j(0) ^ w(0) 

- sin JO) cos i<°> cos ^ ♦ J««» „ ( J sin 2 P) m JO) ^ JO) _ e (0) JO) ^.(0)^ JO) + , 0 e <0) ^2 *(0) ^ ^>^2 w (0) ^ i( 0) ^ JO) m j 0) 

- f * (0 > sin 2 i*°> sir J JJ’ cos i<°> cos f J* cos JO) + 10 e (0) ^2 JO) ^2 J J> sin JO) ^ JO) _ , 5 *<0) si „2 j<0) sin 2 „<£> sin 3 ^(0) cw *(0) 

- 10 P) sin 2 P) sin 3 u'1' sin 2 J°) cos i<°> cos J J* cos c^°> ♦ £ *«» sin 2 i<°> sin 3 J* cos i <°> cos cos gj<°> - IIP) si „2 j(0» ^4 „<£> sin JO) 
cos - 10 P) ^ j(0) sin'* J J> *n 3 w <0) OT j<0) _ «<0) ^ J J> ^ j(0) ^ „<«>> JO) «. e ,0> sin 2 JO) ^ j 0) ^ i(0 ) 

- n<°> i (~ Je (0) sin 2 i<°> cos i<°> cos ♦ P> cos cos w<°>) 4 n<°> t Q (-A P) sin 2 i<°> cos i«» cos co<°> + P) cos *«» cos u><°>) I 

4 ' I 



Element e 


e d) * e ( 0 ) + K 2 (e) [I 2 (e) - S 2 (e) n<°) 1 1 1 + Kg (e) [I 3 (e) - S, (e) n<°) tj_ (4-90) 

l 0 


K 2 (e) = 


K 3 (e) = 


B (0) 4 , 2 

(l-e<°)2)2 

M 

3(l-e<°) 2 ) 3 


4b(°) 6 /i\ . 3 


(4-91) 


(4-92) 


Element iu 


U<1> = “V + K 2 M [I 2 (w) - S 2 (u) n<°> t ) ‘ + K 3 (w) [I, iu) - 8, <«) n<°> tl- (4-93) 

0 <0 


K 2 (<*>) - 


B<0) 4 ,2 


K 3 (u) = - 


e<«) ( 1 — e (0)2,2 

© 


B (0) 6 /ia\ . 3 


^(0) (!— e ( 0)2 } 3 


(4-94) 


(4-95) 


Element i 


i (1) - iq + Ko (i) ;i 2 (i) - S 2 (i) n<°) t]* + K 3 (i) [I 3 (i) - s 3 (i) n (°) t]l (4-96) 

0 t 0 


. B<0) 4 r 2 

K 2 (i) = r e 


k 3 (*)- 


2(i_ 8 (0) 2 )2 

ft)'. 1 


B (0)*>(h\„ 3 


(l- e (0)2)3 


(4-97) 


(4-98) 


Element Q 


0 > + K 2 (n) [I 2 (n)*-S 2 (n)n(°)t]! + K 3 (n) [I 3 (ft)-S 3 (ft)n(°)t]! (4- 

0 t 0 


99) 


K 2 (12) = 


2B<°) 4 r e 2 

7i- e (0)2)2 

B<0) 6 /_ J 3\ . 3 


K 3 (J2) - 


(f) 


(l-e(°)2)3 


(4-100) 


(4-101) 
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Element B 


B<1> “ B< i > + K 2 <») lh <»> ~ S 2 (B) „«» tl‘ + K 3 (B) [I 3 (B) - S3 (B) „«» {] ‘ t 4 ’ 102 ' 


B«» 5 r 2 

K 2 (B) 'e 


,l_e(0)2)3 


(4-103) 


4B<°) 7 /- J l\ 

\ J 2> 


K 3 (B) 


3 (J — e (°) 2 ) 4 


(4-104) 


Element M 


M (1) - M ( K 2 (M) [I 2 (M) - S 2 (M) n(°) t j ‘ + K 


iv + Kg (M) |I 3 (M) -S 3 fM)n<°){]! (4-105) 

0 t 0 


B«» 4 r 2 
K 2 (M) = e 


e (0) (1 _ e (0)2 ) 3/2 


K 3 (M) 


4B(°) 6 /ii\, 3 

— 2 ' 6 

3e (0) ( i_ e (0)2,5/2 


(4-106) 


(4-107) 


As discussed in the beginning of Paragraph 4.3, the second approximations to 
the solutions are then formulated as 


e = + e 

(4-108) 

co = co(0) +e w (l) 

(4-109) 

i =i(0) + ei (l) 

(4-110) 

£2 * £2^ + e $2(1) 

(4-111) 
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B = b(°) + e B^ 1 ) 


(4-112) 


+ (4-113) 

The second approximation for true anomaly V_ is obtained by a Fourier-Bessel 

2 

expansion involving M and e. To the order of e , this expansion is (Reference 6, 

P- 89) 

i'*M + 2esinM +— sin 2 M (4-114) 

4 


where e and M are given by Equations (4-108) and (4-113), respectively. 

As mentioned at the beginning of this section, the super-one solutions due to 
drag are negligible. This fact is illustrated by the following consideration. The 
second approximation to the total solution has the form: 


e-eW + ce* 1 ) 

where represents any orbital element in the set (B, e, i, Q, uu, and M). The eE^ 
terms are short -periodic (see Paragraph 4.3.4) and are composed of integrals of the 
form: 

fi(1) = k 2 ( e ) f [Periodic] d + K 3 (E) f [Periodic] d *^°) + K D (E) / [Periodic] d *X°) 


where K (E) is the constant associated with J effects, K (E) the constant associated 
2 ~2 3 

with Jg effects, and K D ( E ) the constant associated with drag effects. Since 
the above integrands are composed of trigonometric functions which do not 
yield overall solutions secular in the integrated terms will be trigonometric 

functions having amplitudes proportional to the respective constant K (E), K (E) or 

u O 

K d (E). An order-of-magnitude analysis has revealed that K D ( E ) is considerably 

smaller than K (E) and K (E) for each element. Specifically, for a low-eccentricity 
2 3 2 

orbit (e = 0. 0055) and a C^(A/m) of 0. 02 m ,/kg, the relative magnitudes of these 
constants were found to be approximately: 
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E 

£K 2 (E) 

£K 3 (E) 

— 

£K d (E) 

B 

1 x l(f 5 

2 x 10" 8 

0.8 x 10" 16 

e 

1 x 10" 3 

1 x 10~ 6 

0.9 x 10" 14 

ti 1 

8. 8 deg 

0.009 deg 

0.9 x 10‘ 10 

M 

9 . 0 deg 

0.009 deg 

— 

0.9 x 10“ 10 


Because eK (E) is 10 to 11 orders of magnitude less than eK (E) and 8 to 9 

LJ £ 

orders of magnitude less than eK (E), it would appear that drag effects can justi- 

O 

fiably be neglected in deriving the super-one solutions. To verify this, a computer 
run was made for the elements B, e, and uj in which even the super-one solutions 
due to T were neglected. As expected, there was very little difference in the super- 

u 

one solutions with and without the effects of J . Consequently, there would be even 

O 

less difference in the super-one solutions with and without the effects of drag. 

4.3.3 Procedure for Obtaining the Third Approximations to the Solutions 

The third approximation of the solution for any element IS has the form: 

E = + e E^ + e 2 E^ 2 ^ 

( 2 ) 

Procedures for obtaining the E solutions have been established when considering 
oblateness only, but as yet have not been executed to the point of completely deter- 
mining the third approximations. These procedures are outlined in this section. 
Included within their development are the steps necessary to obtain expressions 
for C(t) ; a detailed discussion of these steps is given in Appendix G. 
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The procedure begins by further expanding the basic differential equations of 

2 

satellite motion (Equations (2-34) through (2-39)) to the order of f . These equations 
will have the form: 


_dE_= e f (e(°)) + e 2 g (e(°), E^ 1 )) 
dt 


Specifically, for the element £: 


-f-' e f, e<°>, i«», «<»>, ,«») + e 2 g . (B (0)_ B (l) e <0)_ e (l) j(0) i (l) ) w ( 0 )_ „<!)_ „<(»_ ,(1), 


(4-115) 

which is a functional extension of Equation (4-3). The asymptotic series solution for 
i_ has the form: 

i (t) = i<0) (i) + e i«!» (t, t) + e 2 i< 2 > ft. t) (4-116) 


where 


t = t(l+a' 2 e 2 ) 
t = e t 


Differentiating Equation (4-116) with respect to time yields 



3i(l) i3i (0) 

+ e 2 

di(2) + . diW + djd> 

dt at 

at at 


at at at 


(4-117) 


which is an extension of Equation (4-7). Equating coefficients of like powers of £ 
from Equations (4-115) and (4-117) results in three partial differential equations to 
be solved (two of which are the same as before): 

a i<° ) =0 (4-8) 

a t 


4-39 


3iW^i<0> 

(4-9) 

9 t 9 1 1 


3j(2) + - 

at z a t at 1 

(4-118) 


Solving Equations (4-8) and (4-9) produced second approximations in which it was 
necessary to assume the integration constant C(T) to be a true constant (see 
Paragraph 4.3.2). Solution of Equation (4-118) will permit the functional determina- 

"K (2) 

tion of C(t), as well as a partial determination of i . (Thus, the process of obtain- 
ing the higher-order solution E^ serves to complete the E^ solution.) 

Equation (4-118) will now be considered in more detail. Equation (4-8) implies that 
i^ is a function of 7 only, thus: 


a£2) 

at 



(4-119) 


The solution for i^ is composed of terms nonsecular inj^ and a function which depends 
only uponji, i.e. , from Equation (4-28): 


i( 1 ) = f N (E<°)) + C i (t) 


where 


[ N 


K 2 (i) 


I 2 (i) - S 2 (i) n<°) t 


+ K 3 (i) 


I 3 (i)-S 3 (i)n<°)t 


Equation (4-119) can, therefore, be written as 


dt 9t dt 


(4-120) 
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The solution of Equation (4-120) can be obtained by the same procedure used to solve 
Equation (4-9). Thus: 


t + /g i dt + C' 1 <f) (4-121) 

at at 

where C! (t) is a constant of the jt integration associated with the i' ' solution. At 
this point, it is necessary to again apply the first uniformity condition. Terms 
secular in^ are collected in Equation (4-121) and set equal to zero, yielding: 

d C: (t) - - , bt K _ 

L — t = f gidt -/ ildt (4-122) 

d t sec. sec. d t 


This equation allows the functional determination of C(t), and its solution is discussed 
in Appendix G. 

( 2 ) — 

The V solution thus becomes only a function of terms nonsecular in t_ and the 

integration constant, i. e. , 

i(2) = — / dt + J g: d t + C j (t) (4-123) 

N.S. dt N.S. 


These nonsecular terms can be evaluated by the same technique used to arrive at 

Equation (4-28). Again, the integration constant, C! (t), must be assumed truly con- 

(3) 1 

stant or determined from the E y ' solution. 

( 2 ) 

The procedure for computing the complete E solutions, as outlined above, 

is relatively straightforward; however, it involves many long expressions and taxes 

even the capabilities of FORMAC. It is uncertain at this point what benefits would 

( 2 ) 

be derived, in relation to the work involved, by obtaining the complete E solutions. 
As will be discussed in Paragraph 4.3.4, element solutions are composed of 
short-periodic, long-periodic and secular components. Secular and long-periodic 
terms are the most important in long-term ephemeris prediction, and these terms 
-ire presently contained entirely within the E^ solutions. The short-periodic terms, 
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on the other hand, are contained entirely within the solutions. It is anticipated 

that the C(tjs will add? secularity to the E^ solutions; however, the E^ solutions 

(3) 

will be purely periodic until the E solutions are evaluated, at least to the point of 
~ ( 2 ) 

determining C’(t). Thus, the E terms will, most likely, only add terms on the 
2 

order of £ to the periodic solutions. Since these indications have not been com- 

( 2 ) 

pletcly verified, it is recommended that the nature of the E solutions be further 
investigated before undertaking the laborious procedure of completely solving for 
them. 

( 2 ) 

Undoubtedly, carrying the E' solution procedure to the point of determining 
C(t) for each element would yield certain benefits. For instance, having a functional 
expression for each C(t) would eliminate one requirement for periodically updating 
the epoch values of the elements and associated parameters. (Since each C(t) is 
presently assumed constant, it is one of the parameters that must be periodically 
updated. ) A second benefit would result from the fact that second-order and 
secular effects could most likely be represented in the overall solutions via C(t). 
(The importance of these effects is discussed in Section 5. ) 

4.3.4 Physical Interpretation of the Solution Components 

The analytical investigation of perturbations! effects on a satellite shows that 
(Reference 3, pp. 361-362): 

1. Certain elements experience secular variations from their epoch 
values, as well as periodic variations about these epoch values 

2. Other elements have only periodic variations. 

Earth oblateness, for example, causes secular variations in the elements 0, 
w and M_, and very small periodic variations in all the elements. Similarly, atmos- 
pheric drag causes secular variations in B (or a), e and M, and very small periodic 
variations in all the elements. 
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Among the periodic variations, a distinction is made between long-periodic 
variations (periodic with respect to w or multiples of £) and short-periodic variations 
(periodic with respect to linear combinations of v and w). To visualize these effects, 
consider Figure 4-1. The superposition of all variations depicted in Figure 4-1 yields 
what is referred to as the osculating element. Consequently, the set of osculating (or 
instantaneous) elements defines the continually changing elliptical orbit. 

Inspection of the solution equations for c^, i^, and 

(disregarding, for the purpose of this discussion, the Keplcrian variation in M^) 
reveals that these solutions are secular with respect to^ and/or periodic with respect 
to uj . Consequently, these solution components represent a superposition of the 
secular and long-periodic variations depicted in Figure 4-1, and, as such, represent 
the mean elements. (A mean element is normally defined as the osculating element 
minus the short-periodic variation; however, as discussed in Appendix I, there are 
other convenient definitions for a mean element.) 

Inspection of the solution equations for e^, i O^, and 

reveals that these solutions are short-periodic. Consequently, these solution components 
represent short-periodic variations of the elements about their mean values. 

In summary, letting E denote any element in the set (B, e, i, fi, m and M), the 
physical interpretation of the solution components is as follows: 


E 


E<°> + CE* 1 * 

short-periodic variation about the mean 
mean value of the element 


I osculating value of the element 
( 2 ) 

(NOTE: Since the E solution components have not been derived, they are not depicted 

( 2 ) 

above. As discussed in Paragraph 4.3.3, it is anticipated that the E component will 
add secularity (with respect to the slow time-variable t) to the component through 
the constant of integration C(t).) 
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Figure 4-1. Typical Orbital Element Variations 
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4.4 SYNOPSIS OF THE FORMAC PROGRAM USED rN OBTAINING THE ASYMPTOTIC 
SERIES SOLUTIONS 

As was indicated in Paragraph 4.1, general perturbation methods require a 
great amount of analytical labor in formulating and integrating the equations of 
motion. Equations such as (4-18) and (4-58), for example, involve a large number 
of integrals of the form: 


/ sin p x cos ^ x d x (P, Q * 0, 1, 2, . . .) 


While these integrals are basic, each generally requires several tedious 
recursions in its analytical evaluation. 

In addition, more complex integral forms may arise, such as 

K/ 8 in P xcosQx dx ( n * 1, 2, . . .) 

(1+e cos x)^ 

To alleviate the analytical labor required in performing numerous evaluations 
of both integral forms, an IBM 7094 FORMAC program (IDIGTE) was developed which 
provides the required expansion and integration capabilities. This program consists 
of a FORMAC driver and a set of subroutines which effect the required integrations. 
The driver performs all required manipulations of each input integrand, determines 
the integration parameters P, Q, N and the ’’constant" £ and then transmits these 
quantities to the driver routine of the integration package (the set of routines which 
perform the required integrations). The integration package driver then identifies 
the integrand involved, makes any necessary variable transformations, and calls 
upon the proper subroutine to carry out the integration. 

The complete solution of an integrand usually requires solving several s*jb- 
integrals (special cases), and each integration package subroutine is designed 
to integrate a given type of subintegral. Basically, each integration is carried 
out by substituting the prederived and precoded solution for that particular integral 
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(i. e. , the Integral determined by the values of I> £ and N). These "Integrated" 
results are then transmitted back to the integration package driver (where inverse 
transforms are performed, if necessary), and the results passed on to the FORMAC 
driver for simplification and output. 

A detailed description of this program is provided in Reference 17. 

4. 5 SYNOPSIS OF THE GENERAL PERTURBATION FORTRAN PROGRAM FOR 
NUMERICALLY EVALUATING THE ASYMPTOTIC SERIES SOLUTIONS 

A FORTRAN double -precision computer program (GENPUR) has been developed j 
for the UNIVAC 1108 to numerically evaluate the analytical solutions derived 
in Paragraph 4.3. Currently, the program reflects only those perturbations due 
to earth oblateness (J and_J ) and tangential atmospheric drag, but it is structured 

“6 u 

to readily accommodate additional perturbations, such as higher-order harmonics, 
low-level thrusting, solar radiation, etc. At user option, asymptotic series 
solutions, through the second approximation, can be eva 'ated when considering 
either earth oblateness or the combined effects of oblateness and drag. 

As indicated in Paragraph 4.3, certain assumptions were made in order to effect 
the integrations involved and to evaluate the corresponding integration constants. 

These assumptions appear to be physically reasonable, if they are considered to 
hold over time intervals which are not extreme. With this in mind, the program is 
structured to make use of an updating scheme, whereby the solutions evaluated over 
a gWen time interval (At) are expressed in terms of constants and epoch values of 
the elements (both osculating and mean) computed at the beginning of that time 
interval. These solutions are then used to recompute the constants and epoch 
values prior to the solution evaluation over the next time interval. Included in this 
scheme is a procedure for updating the Fourier coefficients appearing in the series 
approximation to the atmospheric density function (see Paragraph 4. 3. 1. 2). At the 
beginning of each time interval, these coefficients are evaluated by using the 1970 
Jacchia atmospheric density model. 
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Even though the asymptotic solutions do not yet include and second-order 
perturbations, estimates of these effects for the elements w , Q and M were 
temporarily implemented using Brouwer's solutions (see Section 5) to make meaning- 
ful comparisons with actual satellite data. 


A complete description of this program (referred to as the GENPUR program) is 
provided in Reference 18. 
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SECTION 5 - COMPARISON OF RESULTS FROM THE GENERA L 
PERTURBATION PROGRAM 


The asymptotic series solutions through the second approximations have been 
implemented into a computer program (GENPUR). Output of the GENPUR program 
consists of mean elements (E^) and osculating elements (E^ + e E^). (As dis- 
cussed in Paragraph 4. 3. 4 and shown in Appendix H, the mean elements contain the 
very important long-periodic and secular effects, while the osculating elements result 
from adding short-periodic effects to the mean elements. ) The purpose of this sec- 
tion is to thoroughly discuss and compare the results obtainable from GENPUR, 

The comparison of GENPUR results is conducted in two parts. First, the 
validity of the mean element solutions over long time periods is established by a 
comparison with mean elements derivea 'rom Smithsonian tracking data, along with 
corresponding solutions from the MSFC Orbit Lifetime Program. Numerical integra- 
tion programs such as COWELL and ENCKE would have been ideal for this comparison, 
but they are restricted in their application to relatively short time intervals (20 days). 
On the other hand, use of elements derived from tracking data provides the opportunity 
of observing the actual behavior of an orbit, since there are usually small forces in an 
actual environment that are never modeled. In the second part of the comparison, 
osculating element solutions from GENPUR are compared to the results of two numeri- 
cal integration programs, COWELL and SPERTB. These solutions are analyzed for 
only 8 days, since they are merely short periodic additions to the mean elements. 

5. 1 COMPARISON WITH SMITHSONIAN TRACKING DATA 

The Bakci. -Nunn system operated by the Smithsonian Astrophysical Observatory 
(SAO) is a source of very accurate satellite tracking data. The purpose of this 
section is to present a detailed comparison of GENPUR results with SAO tracking data 
for three satellites, namely: Explorer 7, Explorer 1 and SA-5. 
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The Baker-Nunn camera is an instrument with very high accuracy individual 
measurements. The timing accuracy of observations is approximately 0.001 second 
(corresponding to an in-track error of 10 meters for a satellite at a 1000-km altitude). 
Average positions are accurate to within 3 to 4 seconds of arc. The camera takes a 
time exposure of a satellite which is in sunlight, while the camera is in darkness. The 
exposure is interrupted by a rapid operation of the shutter so that the photograph 
appears as a dashed streak of light. The time of the middle interruption is recorded 
with an atomic clock. Appearing on the photograph with the dashed streak (which is 
the satellite) will be point sources of light, which are known stars. The locations 
(right ascensions and declinations) of these stars are accurately predetermined so that 
the photograph provides a recorded history of where the satellite was in relation to 
known references. The processing of these pictures is done with extreme care, requir- 
ing as long as several weeks to get the final results. 

A series of these measurements are then analyzed by the SAO Differential Orbit 
Improvement program (DOI). The DOI program determines, through a least-squares 
procedure, the set of orbit elements that most accurately represents the satellite 
motion during the period of observation. These elements are published for some satel- 
lites in SAO Special Reports. An example is given in Figure 5-1 (taken from Reference 
19) which shows the elements for the initial history of Satellite 1964-5A (SA-5). The 
elements are given in 1-day increments of Modified Julian Date (MJD); however, they 
are not exactly in the form desired. The analysis of this report uses semimajor axis 
a, eccentricity e , inclination^ , right ascension of ascending node jfl, argument of 
perigee cc , and mean anomaly M. Columns 2, 3, 4, and 6 of Figure 5-1 give a,', f2_ , 
and j. in degrees, and M in revolutions. Column 5 presents the history of eccentricity, 
and Column 9 the history of perigee radius in megameters. Semimajor axis is 
obtained from 

r 



Anomalistic mean motion (in revolutions per day) and its first derivative are given in 
Columns 7 and 8. Information pertaining to the number of observations on which each 
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Figure 5-1. Example of Orbit Elements from SAO Reports 


SAO mm element* -- Satellite 196* 5A 1-29 February 1964 



set of elements is based and accuracy of these observations is given in Columns 10, 

11 and 12. 

The first approach taken in presenting the comparison of computed orbit elements 
to SAO elements was to simply plot the histories of SAO elements with a solid line and 
the computed elements with points. In this manner the actual behavior of the element 
could be observed, as well as how closely the mathematical simulation duplicated it. 
However, the angular elements fJ, cc, and M revolve through several hundred degrees 
(and in the case of M, several thousand revolutions). Thus, a small deviation of the 
computed from SAO would be unnoticed, so a second method of presentation is used - 
a plot of the difference of computed minus SAO. These are much more revealing for 
the three angles £), tc, and M, and are the only ones presented for them. 

Also, the computed results of the MSFC Orbit Lifetime Program (Reference 20) 
are shown in the comparisons. This program is indicative of the current state-of-the- 
art in long-term ephemeris prediction, and, as such, provides a standard basis for 
evaluating the GENPUR results. 

As discussed in Appendix I, there are various ways of defining a mean element. 
The GENPUR definition is essentially osculating minus short-periodic, whereas the 
Orbit Lifetime Program and the SAO DOI program use Kozai's mean elements. The 
essential difference is that in defining mean a, Kozai subtracts an additional term (see 
Appendix I). In the following comparisons, this term is added back to the Lifetime 
Program solution for a and to the SAO definition of a so that all are equivalent. 

The GENPUR program is in a developmental state; consequently, it presently 

lacks, among other things, representation of second-order <L and J effects. These 

—2 — 4 

effects are very important for the elements O and a; , and have a slight effect on M. 

To illustrate their importance, orbits of two satellites were simulated by GENPUR, with 
and without an approximate solution for these effects. The approximate solutions were 
obtained from Brouwer's theory (Reference 21) and are in the form of corrections 
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added to the GENPUR mean elements at each time point. These corrections are 
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Figures 5-2 through 5-5 show the GENPUR errors in O and cc for the Explorer 7 and 
Explorer 1 satellites with and without these approximations. (There was little notice- 
able difference in M.) For Explorer 7, the approximations unfortunately increase the 
error in Q from 0. 3° to -0. 7°. However, they decrease the error in u? from a secular 
1. 1° to a random +0.2°. For Explorer 1, the effects are much more drastic. The 
approximations reduce the error in fi from 6. 5° to 0.45°, and in w from -8. 3° to -0.2°. 
Because these second-order _J and_J effects are so important, the Brouwer approxi- 

£ T 

mations given above presently remain in the GENPUR program and will be included in 
all subsequent comparisons. 
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Figure 5-2. GENPUR Error in Ascending Node for Explorer 7 
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Figure 5-3. GENPUR Error in Argument of Perigee for Explorer 7 
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Figure 5-4. GENPUR Error in Ascending Node for Explorer 1 
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Figure 5-5 
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. GENPUR Error in Argument of Perigee for explorer 1 
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5.1.1 Explorer 7 Comparison 


A 344 -day history of the orbit elements for Explorer 7 beginning on 31 March 

1962 was computed by the GENPUR and Orbit Lifetime Programs. Initial mean elements 

and ballistic coefficients for each program are given below in Table 5-1. The m/C^A 

values were adjusted in each program to yield the best overall simulation of the decay; 

2 

in this case, the resultant values were the same (40 kg/m ). (Note the difference of 
1. 1 km in initial semimajor axis due to the definition of Kozai's mean elements used by 
the Lifetime Program. ) 

Table 5-1. Initial Conditions for Explorer 7 


MEAN ELEMENTS 

GENPUR 

ORBIT LIFETIME 

a 

(km) 

7193.0 

7191.9 

e 


0.03545 

0.03545 

i 

(degrees) 

50.305 

50.305 

O 

(degrees) 

344.40 

344.40 

tc 

(degrees) 

232.44 

232.44 

M 

(degrees) 

179.46 

179.46 

m/C D A(kg/m 2 ) 

40.0 

40.0 


The histories of semimajor axis and eccentricity are shown in Figure 5-6. The solid 
line is a connection of each SAO element point (given at 4-day intervals). The 
asterisks represent simulation results from the GENPUR program, and the circles 
are results from the MSFC Orbit Lifetime Program. Both simulations are nearly 
coincident for a_ and e., and both show extremely good agreement with the SAO elements. 
(Recall that the output of the Lifetime Program and the SAO values of semimajor axis 
nave been adjusted to remove the Kozai correction. ) Semimajor axis decays only 
slightly (0. 5 km) during this interval, so that the orbit is essentially free of significant 
drag effects. The long-period variation in eccentricity due to J is very evident, having 

u 

a period of approximately 110 days and an amplitude of 0. 0008. Note also the relatively 
rough nature of the tracking data, especially for £. There seems to be bad tracking 
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Figure 5-6. Semimajor Axis and Eccentricity for Explorer 7 
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data points in the values of a_ at 48 and 252 days. (These characteristics are evident 
for each element and each satellite. ) 

Figure 5-7 shows two different types of plots. The top figure presents SAO 
tracking values of inclination, along with individual values of inclination from the 
GENPUR and Orbit Lifetime Programs. Tracking values of inclination are fairly 
rough, at least on the scale being usee. The Lifetime Program holds inclination con- 
stant at the initial value (50. 305°). The GENPUR program simulates the secular 
change in inclination due to drag, not evident in this figure, and the periodic change 
due to J which can be seen. (The advantage of having inclination vary is not apparent 

tJ 

for Explorer 7, but will be for Explorer 1. ) 

The bottom half of Figure 5-7 shows the error in ascending node produced by 

each program, i. e. , computed value minus SAO value. Again the results of each 

program are nearly identical and both show fair agreement with the SAO elements. 

There is a secular buildup of error in ascending node to -0. 7° for each program. 

(Recall that the GENPUR program uses Brouwer’s equations to approximate the J and 

second-order _J effects in O, cr, and M. ) 
z 

In-track position of a satellite is primarily a function of argument of perigee 
(U!) and mean anomaly (M). The mean anomaly typically undergoes 5000 revolutions 
in 340 days. It is extremely sensitive to small changes in semimajor axis. For 
example, an error of only 0.4 km in semimajor axis can result in an error of 170° 
in mean anomaly after 340 days. Mean anomaly is very sensitive to gravity and drag 
perturbations; thus, it provides a significant measurement of the accuracy of a 
simulation. Figure 5-8 shows the errors of both programs in cc and M. Both have 
almost identical simulations of cc with no apparent secular error, but only a random 
error of +0.2°. These differences may, in fact, be due to limitations on the accuracy 
of the tracking data, rather than inaccuracy in the simulations. The simulations of 
mean anomaly are somewhat different. The Orbit Lifetime Program shows a periodic 
and secular error buildup of nearly 75°. The GENPUR program, on the other hand, 
exhibits only a secular error buildup of 40°. 
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5.1.2 Explorer 1 Comparison 


A 356 -day history of the orbit elements for Explorer 1 beginning on 2 January 
1964 is shown in Figures 5-9 through 5-11. Elements derived from SAO tracking data 
are shown along with computed solutions from the GENPUR and Orbit Lifetime pro- 
grams. Initial mean elements and ballistic coefficients for each program are given 
below in Table 5-2. Again, the m/C^A values were adjusted in each program to yield 
the best overall simulation of the decay; in this case, the resultant values were 
different . It is thought, that the reason is due to the fact that short-periodic pertur- 
bations in altitude were not considered when determining the Fourier coefficients of the 
GENPUR program. (The difference between the two definitions of a amounts to 5. 0 
km for this orbit. ) 


Table 5-2. Initial Conditions for Explorer 1 


MEAN ELEMENTS 

GENPUR 

ORBIT LIFETIME 

a 

(km) 

7368. 14 

7363. 14 

e 


0.08747 

0.08747 

i 

(degrees) 

33.198 

33.198 

Cl 

(degrees) 

34.01 

34.01 

LC 

(degrees) 

151.27 

151.27 

M 

(degrees) 

50. 112 

50.112 

m /C D A (kg/m 2 ) 

22.28 

25.0 


The histories of semimajor axis and eccentricity are shown in Figure 5-9. The same 
plotting symbols as before are used, so that the straight line is a connection of SAO 
elements, asterisks represent GENPUR results and circles are Orbit Lifetime results. 
Both simulations are nearly coincident for a. and show reasonably good agreement with 
the SAO elements. The simulations are initially about 0.8 km higher than actual and 
then fall about 0. 6 km below actual after 350 days. The reason for this behavior is 
due to omission of daily values of solar flux F^ ^ and geomagnetic index A^ in the 
density model. (The Lifetime Program, when using the daily values, showed nearly 
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Figure 5-11. Errors in Argument of Perigee and Mean Anomaly for Explorer 1 




perfect agreement. ) As yet, input of daily ^ and values is not available for 

GENPUR; hence, both programs were run in a simulated preflight condition using only 

mean values of F,* _ and regression values of A . Note that this orbit is affected 
10.7 p 

considerably more by drag due to the lower perigee than was the orbit of Explorer 7. 
Semimajor axis decayed 24 km, rather than the 0. 5 km for Explorer 7. 

The lower half of Figure 5-9 depicts eccentricity. Both simulations agree well 
with SAO elements, but they are not coincident. Long-period effects of_J are again 

O 

clearly evident with a period of 48 days and an amplitude of 0. 006. A secular decrease 
in the magnitude of e^due to drag *s also noticeable. 

The upper half of Figure 5-10 shows the computed simulations and SAO values of 
inclination. The long-period variation due to J is clear, and is reasonably well simu- 

O 

lated by GENPUR. In the Lifetime Program, however, inclination is held constant. 
Therefore, the GENPUR program shows a significant advantage over the Lifetime 
Program in simulating inclination. 

The lower half of Figure 5-iO shows the errors of GENPUR and Orbit Lifetime 
in simulating ascending node. Both programs show very similar results, having 
maximum errors of 0. 45°. 

Errors in the critical in-track angles a; and M are shown in Figure 5-11. 

GENPUR results are better than the Lifetime Program for <£. GENPUR errors grow 

to a maximum of only 0. 3° whereas Orbit Lifetime errors in cc grow to 0. 6°. Errors 

in mean anomaly for GENPUR are smoothly varying with a maximum of -75°. Maximum 

error in the Lifetime Program is also -75°, but note the peculiar periodic nature that 

it exhibits (which was also evident in Explorer 7). The error in mean anomaly from 

the GENPUR program is easily explained in terms of the error in semimajor axis. 

Simulations of mean anomaly are very dependent upon an accurate value of a. (Recall 

that there were small errors in the simulations of a for Explorer 1, Figure 5-9.) 

Initially, computed _a was too large, which means theoretically that the orbital rfiean 
3 1/2 

motion (equal to (ji/a ) ) would be too slow and mean anomaly would not change as 
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rapidly as it should. Figure 5-11 shows that this is, in fact, what actually happened. 
The GENPUR value of mean anomaly initially falls below actual. Then, as the computed 
a becomes close to the actual a at 260 days and falls below actual at 275 d?.ys, the 
error in mean anomaly levels off at the maximum -75° and returns to only -18°. There- 
fore, had the GENPUR simulation of a been better, the error in mean anomaly would 
have been much less. 

5.1.3 SA-5 Comparison 

A 334 -day history of the orbit elements for the SA-5 satellite beginning on 1 
February 1964 is shown in Figures 5-12 through 5-14. Elements derived from SAO 
tracking data, at 2-day intervals rather than the 4-day intervals of the previous satel- 
lites, are shown along with computed solutions from the GENPUR and Orbit Lifetime 
Programs. Initial mean elements and ballistic coefficients for each program are given 
in Table 5-3. (The difference between the two definitions of a amounts to 5. 64 km for 
this orbit. ) 

Table 5-3. Initial Conditions for SA-5 


MEAN ELEMENTS 

GENPUR 

ORBIT LIFETIME 

a (km) 

6889.68 

6884.04 

e 

0.0358 

0.0358 

i (degrees) 

31.4561 

31.4561 

Q (degrees) 

161.797 

161.797 

co (degrees) 

150.01 

150.01 

M (degrees) 

34.56 

34.56 

m/CpA (kg/m 2 ) 

88.22 

106.0 


The histories of semimajor axis and eccentricity are shown in Figure 5-12. 

Plotting symbols and notation are the same as before. Both simulations are nearly 
coincident for the element a, but neither agrees very well with the SAO values. The 
simulations agree reasonably well for the fiist 60 days, but rise above actual by 2 km 
at 100 days and then fall below actual by -3 km at 334 days. This error was encountered 
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Figure 5-13. Inclination and Errors in Ascending Node for SA-5 
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in a previous study and is believed io be due to an inaccuracy of the 1970 Jacchia 
density model at lower altitudes. Since the semimajor axis is not well simulated, it 
is anticipated at this point that there will be relatively large errors in the simulations 
of O, tc, and M. Of the three satellites investigated, drag had the most significant 
effect on SA-5. The initial perigee altitude was only 270 km (versus 340 Ian for 
Explorer 1 and 560 km for Explorer 7) so that the atmospheric density at perigee was 
significantly greater than for the other satellites. In 334 days, the semimajor axis of 
SA-5 decayed by 64 km versus 24 km for Explorer 1 and 0.5 km for Explorer 7. 

The lower half of Figure 5-12 depicts eccentricity. The two simulations are 
nearly coincident, but again do not agree well with SAO values. The reason is the 
same as lor the discrepancy in a. Long-period effects in e^ are clearly evident with a 
period of 36 days and an amplitude of 0. 0008. A secular decrease in the magnitude 
of e due to drag is also noticeable. 

The upper half of Figure 5-13 shows computed simulations and SAO values of 
inclination. Long-period variations are not evident in the SAO values. In fact, the 
random fluctuations in the SAO data are larger than the amplitude of the long-periodicity, 
implying that the resolution of the SAO data was not accurate enough to show the long- 
periodicity. The SAO elements also show a very interesting phenomena at 150 days, 
where the average value of inclination seems to change from 31. 456° to 31. 465®. It is 
liard to imagine what physical force could cause this change other than a powered 
plane -change maneuver; however, no such maneuver was performed by SA-5. It can 
be concluded that the GENPUR simulation of inclination for SA-5 is as accurate as the 
SAO elements. 

The lower half of Figure 5-13 shows the errors of GENPUR and Orbit Lifetime 
Programs in simulating ft. The two programs agree with each other for the first 210 
days, but then the errors diverge. The error of the Lifetime Program decreases more 
rapidly than does that of the GENPUR program. The reason is that, at this point, the 
Lifetime Program simulation of a falls slightly below that of GENPUR. Neither pro- 
gram shows particularly good agreement with SAO values, both having a maximum 
error of nearly 0. 7°. This was as expected, since a was not well simulated. 
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Errors in the critical in-track angles u? and M are shown in Figure 5-14. The 
two programs are not coincident in simulating a.', but it would be hard to say which is 
better. Both show errors ranging from -1. 0° to +1. 2°. Again the trouble is due to a 
poor simulation of a. 

Errors in mean anomaly are very dependent upon the simulation of a. Thus, the 
poor simulations of a by both programs are very evident in their large errors in M. 

The GENPUR error in M ranges from -80° to more than +180°. The Lifetime Program 
error ranges from -135° to more than J -180°. Simulations of the elements of SA-5 
clearly demonstrate the importance in orbit ephemeris prediction of having a good 
simulation of semimajor axis (which depends upon the use of an accurate density model). 

5.1.4 Summary of Tracking Data Comparisons 

More than 300 orbit days for Explorer 7, Explorer 1 and SA-5 have been simu- 
lated by the GENPUR and Orbit Lifetime Programs. A summary of the errors of the 
simulations for each orbit element and each satellite is shown in Table 5-4. (Recall 
that the GENPUR program does not yet contain asymptotic expansion solutions for 
second-order secular effects of ^ aid J , but uses Brouwer’s equations. ) A + sign 

Z “r 

indicates that the error was more or less random, and is the type desirable for all the 
errors. A single number means that the error steadily increased to the value given, 
whereas two numbers indicate that the error grew to the first number and then reversed 
direction and attained the level of the second number. For example, the error by 
GENPUR in a for Explorer 1 first grew 0. 8 km and then reversed direction to -0. 6 
km. 

In general, the errors of both programs are nearly equal with one or two excep- 
tions. Having the J effects on inclination included in GENPUR results In only a 

u 

-0.003° error rather than the -0,008° to 0.004° error of the Lifetime Program (for 
Explorer 1). The argument of perigee for Explorer 1 was better simulated by the 
GENPUR program, as was the mean anomaly for Explorer 7. 
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Maximum Error 
During the Full 
Simulation in: 


EXPLORER 7 


(deg) 

(deg) 

(deg) 

(deg) 


GENPUR 

Lifetime 

GENPUR 

±0.3 

±0.3 

0. 8, -0.6 

-0.0002 

-0.0002 

±0.0003 

±0.01 

±0.01 

-0.003 

-0. 7 

-0. 7 

0.45 

±0.2 

±0.2 

0.3 

40 

70 

-75 


EXPLORER 1 


R Lifetime 

6 0 . 8 , - 0.6 
- 0.0002 
-0.008,0.004 
0.45 
0.6 
-75 



GEN PUR 
2. 0,-3. 0 
-0.0006 
- 0.01 
0. 7 

- 1 . 0 , 1.2 

-80, 180 


Table 5-4. Summary of Errors by GENPUR and Lifetime Programs 
in Simulating Orbits of Three Satellites 


Lifetime 
2.0, -3.0 
-0.0006 
- 0.01 
0.7 

- 1 . 0 , 1.2 

-135, 180 


Note: The GENPUR simulations 
contained Brouwer's approxima- 
tions of second -order J and J 
effects. 4 




5.2 COMPARISON WITH NUMERICAL INTEGRATION PROGRAMS 


Once an accurate history of mean elements is available, the osculating elements 
can be obtained by adding short-periodic terms. These short-periodic terms are 
primarily functions of the mean elements. It is not necessary to verify osculating 
elements for long time intervals, providing the mean elements are good. (If the 
osculating elements are good for short periods of time, they will be good throughout 
any given interval providing the mean elements remain satisfactory. ) 

SAO tracking data do not contain osculating elements. Therefore, a different 
method of comparison was necessary to verify osculating element solutions. The 
MSFC COWELL and SPERTB numerical integration programs use osculating elements 
exclusively; therefore, Table 5-5 shows a comparison of GENPUR osculating elements 
to those from the COWELL and SPERTB programs. Initial conditions are the initial 
SAO elements for Explorer 7. A period of 8 days was simulated, so that the GENPUR 
mean elements experienced little error. 

In Table 5-5, the osculating element solutions are shown at the end of 1 day and 
8 days. Four simulations were run, namely: the COWELL program, the SPERTB 
program with and without effects, and the GENPUR program without (and second- 
order^) effects. The GENPUR results are within the differences between the COWELL 
and SPERTB programs; thus, the GENPUR osculating element solutions are excellent. 
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AFTER 1 DAY 

AFTER 8 DAYS 

; COWELL 

SPERTB 

SPERTB 
(w/o J 4 ) 

GEN PUR 
flv/o J 4 & j|) 

COWELL 

SPERTB 

SPERTB 
(w/o J 4 ) 

r 

GENPUR 
(w/o J 4 & £) 

0.035699 

0.035698 

0.035698 

0.035679 

0.035456 

0.035460 

0.035453 

0.035432 

235. 734 

235.733 

235.733 

235.713 

260.436 

260.352 

260.428 

260.353 

50.299 

50.299 

50.299 

50.299 

50.304 

50.304 

50.304 

50.305 

340. 188 

340. 183 

340.188 

340. 192 

310.865 

310.869 

310.865 

310.900 

7192.741 

7192.735 

7192. 737 

7192.812 

7194.044 

7193.997 

7193.980 

7194.032 

263.862 

263.861 

263.866 

263.827 

133.290 

133. 4 oO 

133.459 



133.377 


Table 5-5. Short Term Comparisons of Osc ulating Elements 













SECTION 6 - CONCLUSIONS AND RECOMMENDATIONS 


The basic objective of this research project has been to develop, through applied 
research in general perturbation theory, perturbation techniques that provide an 
accurate and rapid long-term ephemeris prediction capability for satellites in earth 
orbit. The approach taken was to use two-variable asymptotic series in obtaining 
approximate solutions to the Lagrange planetary equations of orbit motion. This 
technique constitutes a relatively new approach to the ephemeris prediction problem 
and, while it is not yet on a rigorous mathematical basis, offers several potential 
advantages (as discussed in Paragraph 4. 2). In this study, it was found that two- 
variable asymptotic series can be successfully applied to the problem of artificial 
satellite motion under the combined influence of gravity and drag. The first and 
second approximations of element solutions derived by asymptotic series agree in 
form to those derived in other established theories. 

Of the potential advantages which the asymptotic series method offers, two were 
found to be of significant aid thus far. Since the method employs two time scales, 
the solutions obtained tend to group naturally by physical effects, i. e. , they group 
into secular, long-periodic and short-periodic components. Therefore, it is not 
necessary to use a procedure such as Kozai’s in which the disturbing function is 
resolved into secular, long-periodic, and short-periodic parts. Second, the error 
involved in a given series approximation is of the order of the first neglected term. 
Consequently, the asymptotic solutions are naturally structured to include the 
dominating effects of each perturbation in the initial approximation. Furthermore, 
a control of the expected error is provided by selection of the expansion parameter e. 

Currently, the asymptotic eolutions have been obtained through the second 
approximations when considering earth oblateness and tangential atmospheric drag. 

In these approximations, it was found that the E^ solutions contain the very important 
secular and long-periodic effects, while the E^ solutions contain the short -periodic 
effects. The E^ solutions were derived by first obtaining simultaneous solutions to 




the differential equations for the elements e^ 0> and these solutions were then used 

to obtain solutions for the remaining elements. These solutions were carried only 
through the first power of eccentricity. Because of the importance of the com- 
ponents m the total solution, it is recommended that they be investigated further. 
Specifically, extension of the simultaneous solutions to include more elements and 
retention of higher orders of eccentricity are recommended. 

As indicated in Paragraph 4. 3.2, functional forms of the integration "constants” 
C(t) have not been analytically determined, and are currently evaluated by use of an 
updating procedure. If the functional form of each C(t) was available, one require- 
ment for using the update procedure would be eliminated; furthermore, the second- 
order secular effects of oblateness perturbations would be contained in these functions 

(see Paragraph 4.3.3). Analytic determination of these "constants" requires partial 
( 2 ) ( 2 ) 

development of the E' ’ solutions. It is anticipated that the E v 1 solutions, themselves, 

will be purely periodic except for their integration "constant", C'(t). Therefore, it 

( 2 ) 

is recommended that the E solutions be investigated, at least to the point of deter- 
mining the functional forms of C(t). 

During the study, it was found that some form of automated manipulation capa- 
bility is absolutely essential to the accurate and timely solutions of the equations 
involved. Many operations on very lengthy expressions are required, such as expan- 
sions, integrations, substitutions, simplifications, etc. Furthermore, an automated 
method for uniform presentation of results is highly desirable. Therefore, the 
FORMAC language was used to write a computer program that performs these opera- 
tions and presents the results in a convenient manner. As a result, a great deal of 
experience was gained in the use of FORMAC; and limitations of the language, such as 
lack of identity recognition, core storage requirement, problems in subroutine com- 
munication, etc. , were encountered. (A thorough discussion of these problems is 
given in Reference 17.) The necessity of this automated manipulation capability in 
providing accurate and timely analytical results cannot be overemphasized, and the 
development of the FORMAC program is considered to be a major accomplishment 
of the project. 
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The second approximation solutions of orbital motion using, two -variable 

asymptotic expansions have been implemented into a UNI VAC 1108 computer program 

(GENPUR). A comparative study of the results obtained using this program showed 

it to be very accurate, especially when Brouwer’s approximations (Reference 21) of 

the second-order J and J, effects are used. (As yet, these effects have not been 
-2-4 

determined by the methods of asymptotic expansion.) For example, errors in the 

solutions for the short-period effects in mean anomaly for Explorer 7 were less than 

0. 09 degree during 8 days, which was less than the difference between the standard 

COWELL and SPERTB (Reference 10) numerical integration programs. Furthermore, 

errors in the long-term solutions (i. e. , mean element solutions) were generally less 

than or equal to the errors of the MSFC Lifetime Program (Reference 20). (These 

errors in the long-term solutions may possibly be reduced when the approximations 
2 

of ana J ^ effects are replaced by the asymptotic series solutions. ) 

The run time required for an ephemeris prediction program is always of utmost 
importance. The GENPUR program is extremely fast and has the potential of being 
even faster. For example, the run time required for the simulation of the Explorer 7 
satellite over a 360-day period was 94. 5 seconds when using an update interval of 24 
hours. Increasing this interval to 96 hours resulted in no noticeable loss of accuracy, 
and the run time was reduced to only 23.4 seconds. In comparison, the run time 
required for the same orbit using the MSFC Lifetime Program (with a 2-day step) 
was 148. 9 seconds. 

Even in its present developmental state, the GENPUR program has clearly 
demonstrated the soundness of the approach taken herein to compute long-term satel- 
lite ephemeris. Before being placed in a production status, however, there are 

certain additions to the program which should be made. Errors in the element solu- 

2 

tions for w, and M could possibly be reduced by an accurate representation of J 
and effects. Even with the approximations now being used, the GENPUR error is 
comparable to, or less than, that of the Lifetime Program. Since two-variable 
asymptotic series represents a different approach to ephemeris prediction, it is 
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quite possible that this method could result in more accuracy than existing solution 
methods. 

Another addition recommended for GENPTJR, which could result in much faster 
run times than even the 23. 4 seconds mentioned previously, is the use of analytical 
expressions for the Fourier coefficients. One innovative feature of the GENPUR 
technique that contributes to its speed has been the use of Fourier series expansions 
to represent drag effects. The Fourier coefficients are presently determined by use 
of the 1970 Jacchia density model at frequent intervals. If the variations of these 
coefficients for periods of 20 or 30 days could be established analytically, a run time 
of only 6 seconds would be a possibility. Furthermore, the successful development 
of such a model would represent a significant advancement in the state-of-the-art of 
satellite ephemeris prediction. 

An increase in the flexibility of the GENPUR program is also recommended. A 
wide variety of input coordinate systems, as provided in the Lifetime Program, 
would be advantageous. The satellite physical characteristics (mass, drag coefficient, 
and area) must now be held constant in the program. Providing input options for these 
items which allow variations with time and/or orbital position would be extremely 
useful in orbit analyses. Also, it would be desirable to have an input option for daily 
values of solar flux and heating parameters. This flexibility could be easily achieved 
within GENPUR by incorporating many of the corresponding routines of the Lifetime 
Program. 

Once the GENPUR program has been extended as recommended above, it will 
represent an even more valuable tool for conducting astrodynamic investigations. For 
example, King-Hele has stated that the upper atmosphere rotates at a faster rate 
than the earth, but other investigators have failed to confirm this finding. By using 
the GENPUR program to study the long-term evolution of inclination for various orbits, 
an independent estimate of upper atmosphere rotation could be made. Another problem 
which has received little attention is the exact nature of the final decay of eccentricity. 
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It is well-known that an eccentric orbit becomes nearly circular before its ultimate 
decay, but whether it becomes zero or reaches a limiting value is uncertain. Further- 
more, because of its extremely fast run time, GENFUR is ideal for parametric 
studies to identify characteristics of various classes of orbits to aid in mission 
planning activities. 

In summary, this study has demonstrated the successful application of two- 
variable asymptotic expansions and the automated manipulation capabilities of 
FORMAC to the satellite motion problem. The resulting GENPUR computer program, 
although in a developmental state, has clearly exhibited the potential of being more 
accurate and much faster than any existing long-term ephemeris prediction program. 
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APPENDIX A - DERIVATION OF THE 
PERTURBATIVE VARIATION EQUATIONS 


To illustrate the procedure for obtaining Equations (2-7) through (2-12) by the 
method of perturbative differentiation, the equation for e will be derived. 

The polar equation for an ellipse is 


r = 




M- 


(l +Cc+->'J) (|+£umv) 

where the specific angular momentum is given by 


(A-l) 


k 


2 ' 


r v 


(A-2) 


Taking the dot-derivative of Equation (A-l) results in 

r - ^ e. v -^v.-y 
a (i+ e c^.v) a 

or, after substituting Equations (A-l) and (A-2) 


r = ^ e (A -3) 

b 

Substituting Equation (A-2) into Equation (A-l) yields 

3*2 

V V I 

e cev,^ - — — - i 

It is now necessary to take the grave-derivative of this expression, 
remembering that r = 0 (see Paragraph 2.2). Thus, 


e - v v - — 




Similarly, substituting Equation (A-2) into Equation (A-3) yields 


(A -4) 


which becomes, after taking the grave -derivative 

+ 4 r v)- r f y e + LiL /IS (A-5) 

/ Alp vr/ 

Multiplying Equation (A -4) by cos V and Equation (A-5) by sin V and then adding the 
results yields 

F’ ' v^f v 




Since 


e c^-jv 




this becomes 


e ' = ^ FT r(-r + Fvte] 

(NOTE: Equation (2-8) agrees with Reference 2, p. 247) 
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The perturbative variation equations for the remaining elements can be 


derived in a similar manner. 


APPENDIX B - ASYMPTOTIC EXPANSION OF TRIGONOMETRIC 
AND EXPONENTIAL FUNCTIONS 

To illustrate the procedure for asymptotically expanding trigonometric and 
exponential functions, the functions 

0-*X X (I- cT* 

x (l+e.c^K) 

will be considered, where x is used to denote any angle element. 


Begin by assuming the asymptotic series 

X - 6-X^-f 6**X (S + ... * X^V b 


b 


COS X 




But 


Thus, 


=• I--5- £*!&)*= I~*z + ... 

fe a x®+... 

= 6m, X®* fi^X^X®] + 6*f ]* — 


sin x 


= e^x^x^x*!) 

* + 7+... 
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. 2 
sin x 


Thus, 


Note that sin x has the form 

X “ CVj 4 0^6+ cv 3^ 

X — 4 Q.J t J * ^G-i + tX^fc^ -+ ^^a.|4CX 3 6^^56 +G.j£ 

= a., 2 + 6^, *z] + 6*C ]+-•• 

6 ] +6: Z C J +. .. 

= AJ^\ C ''+6fi®4U1?*&]+t X { J+... 


, i 3 

(1 4 e cos x) 

From tiie previous expansion for cos x, 

(I-K2.U ^ x) ” l-h (e®+ £ ^■4...)(6^X <S) 4- 6 /P X^^ X»]) 

- I -#■ e^^x®-4 

= [l+ e^c+ot^] + e[e®Or?x®- e^x^x^J J-f. 


Note that (1 + e cos x) has the form 

^ I 4 £ C*-t 7 x) * C^l 4- O. ^ 6 4 ^3 4r 


fl4 eo^x) ? -^fo.,4 a a e)^a,6 l 7 ts (o> t + & 0 *!f + 5(a,4a 5 fc) r a. 3 
~ ^t^4 £ pkcx," 7 ' CX 2 ~j 4- £ £ P] 4 , . . 
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Thus, 


(l + e^*y = ( l+ x f ’‘ ] % e(xu A,#) Ye ^ A 

+ «Y ] + ... 


d-eV 7/2 

= I- fe A % 2>e r °'A + tV* 1 )- sfe^-te^fcY 9 - t'/ )2 

= (/-/"” J + tf-Se&e 0 '] +&'£'*>'- 2e ( V^... 

Then, 

#/-«»*} + + tY- 2 e ft c S 'J}]' lA 

= fi-ito'y % +£Q(l-t?'y* Ytfs^»e®J+ 6 t ^«®- 2 ^ 1 e« k <J) 

= p./^rK 6 t >... 
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APPENDIX C - PRESENTATION OF THE INTEGRATION RESULTS OBTAINED 
FROM THE FORMA C PROGRAM WHEN CONSIDERING EARTH 

OBLATENESS (J AND J ) AND DRAG 

2 J 


As mentioned in Paragraph 4. 3. 1, a FORMA C computer program is used to solve 

integrals of the form given by Equations (4-18), (4-19), (4-52), (4-53), and (4-58) 

for each element. The program prints the total integrated results L, I , and I_ 

( 0 ) 2 ” 3 U 

for each element, as well as the secular (with respect to v ) parts S^ SL, and S n 
(see Equations (4-20), (4-21), (4-59), (4-60), and (4-61)). These answer arrays are 
presented on the following pages. As discussed in Paragraph 4. 3.2, it is not neces- 
sary to consider drag when formulating the second approximations to the solutions 
since the super-one solutions due to drag are negligible. Consequently, I_ D for eacl1 
element is not required and, therefore, is not presented. 

To maintain consistency with the assumptions made in solving the set of ordinary 

differential equations having t_ as the independent variable (see Appendices E and F), 

the arrays generally include only terms through the order of e (i. e. , terms on the 
2 

order of e , or smaller, are ignored). 

Recall that the numerical subscript on I_ and S indicates the earth harmonic under 
consideration, the subscript D indicates drag, and the parenthetical (x) indicates the 
element x. For example, 

I (0) = Total integrated result for the element f) when considering 

L* 

the second harmonic (J^) 

S (e) = Secular part of the total integrated result for the element e 
when considering the third harmonic (J ) 

J 

S^(e) = Secular part of the total integrated result for the element e 
when considering drag 

Also, recall that in the array for each element, a^ and b^ are the Fourier 
coefficients appearing in the atmospheric density function approximation (see 
Equation (4-57). 
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Element B 


- C J»+ u^-/ v J + Z^'c^v^uPu.-fl*. J 

, . 1 . (?) » ■*• C*J a <5> . a £> 

4^ c V — 1/4^* £ 4^v% V 4^ CO 


4(fi) = o 


^/e) = e JH- liic-V^v^^-L^w 1 

* v r “U w®_ /^^^Vk-«®-* ODA^C&^-JkJ^l 

- ^ ^^L.n.'V W®- <A-.v®+ w® 

j_ l£. .- 3„(S) -** (?) . (5> (?) S ' 1,(o) . 1 A - “S (®) fit) S' *.6>. * (»>- a ,fe> _(o\ 

'~Tf 4*** 4 4*" \) /i«s Co d*o V — a4^*» t 4 ^ v co 4 » ^ V — S/u*~. c V 4*<i w co -* 

, f- - •? . £> - "7 (Sj <o\ i* - -?„£).. 9 (oj ^ . (o) 

+ ^ /ikv, t \) Cer~> Co -+• ^ /4-*-v. (_ aLv. Co 4#^V 


.<•> 


&) = O 

5/S)= 
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Element e 


+ £ ®»C. v W I« J 5 '- 

-u . a.6>, ft , (3 <6j OLa' x '*° - ' \ ft (5) J2-.' 7 .'V. .ft 

7A-*. t 4^—V /U~.c*y C4* to + '4*^- t /4*^v> A»— to c*-> v — J /1a -'- <- <4*~ V 0>v 

.*_ ^ /•j_ jL , ,(•> 


"» -ft. “3 ft . 


r«> 


ft ^/ft , a ft 

3 


+ c 4-C— 4 C—L 0 Co-jto — <§r 4*-~ f 4^v.*(o' £*-» v — .<*»■> i 'c&-^\) '•+ 


^2 (V) - o 


2^fe) = e 


(P>ffj_ „ ft ft IT. , . 

/ ^/d**-t t*-*co — "ZfACvi. 


.ft . i ,0»> „ ft ft r? . ,ft. * ft ft 

# /tu. yP V — iy 4 m» t /4^v V CO 


, ✓ _ (o; , fy L /jtl „ ■?. ft . ft . \ (»\ (•) __ . ?„ft . i (6j - ft. .ft 

4" -ftjX— t >*>/*- Co Co-> V> — 2 ■'**'" < to £#— < to — 2?4L>. c /ti> \) to 6*-* \) 

. /£( . . 7 -ft . 1 ft - ? ft ft 304’ - ?„ft - ,fe) - l (ft (o) i<" - 7.ft „ 7 ,ft. , ft 

-/- ^ -d-*— - £ 4 a*- V -t— Co C<r^\3 *f ty Av*» t V) 4**. CO Cfr-?Co — -ij' 1 4 a*- V ttn to 


4*~ V 

?<_„■ 7-(o) , W fe) 


ft. . jft , *-<•>. Vi ft - ?, .ft. ft A,-, . 7„<^ - < /ft 

j V — 2Ti«< l 4C» \? /Jo*- to V — -? td* a* t .1~~- 'y o 


i (•) 
id**— to C*-r?CO 


ft 


, -2J2. . - "*,(•) „ ^ (ft (»J 7.(o) _ ("i) (ft - . 7.(<0 . T (ft, /*» ”} 

' ty Ao*- j, 4 a*-. v> to — /4**» t /dov-CO (Co-* V — /4 a-~ t xJa— to J 

_ ^ . ,0 - .(4) (•), /£l. (•> . ® . -G>) , £*> C»i (•)$_*- - A > -' l ( *> - f Q ) 

~tj \ 4*-*. t t*— 7 CO -f 75 V '^ Av ' t t^>CO -+■ «y 4**-t /^O— V? C#-^> Ca^> to — Cy ^ A *>-t /1 a— ”V) /J. — to 

(Z' , 7,® . T ft 

” "o” 4ov t A*-V A— to / 


.ft. 






K .<.(').! (*) . '°i 

___ t<r-7 v'^CA-? Co '" -— t ^l-a‘S>'^C 6-*'0 ' Le-->Ui -— /t*-— t 4^- to 

i 3 j 5" ' 7.ft . t ft) t iii 2< , ,7 . r 3 -ft . 7 to) (•) 

-q' A**~ c 4 ^a- \) 4 ^. to -f “j" 4—*- t 4— \) /iA-v- to V’ 6*-3 1* 1 — t /4 a*- V? Cir-sV t<TT> 


3-<'^ - v ft 

/iU*. , t 4**. V 


A—to - 


f'o) ^>*T . 3 .(^O) . ty (i) . * To) 


/4a-w c V 4^* 


^ (e) « C~ Y ^ 60 ^-4- t ‘°^ ) 

4 -|-t> 5 )-f (r-k*»+ i^*)e (oi 
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Element i 


(•) 

4 m- u £*-> V 


^ 6 ) * e 2 £<W^. 

- ^X-3c*4ijJ 4 C'-, 

+ C? /^al— *to’°^ 


4M)~o 


fix <fo 

£*-*CO 


J s/\ ■»•>> .fo to G> .( 

4**- C ft* i ft*co O- \J ft?c 

,_ 7 -ft) ft . 7 fo ,ft fo ft) £ - 7 .^) , .ft . fe> V1 ^ fo) 

“" / CU-. c. /2X— N) ft— CO ft-*t ft* V C#-»CO C 4-vV) £0 » c ft-*V ft? CO 

/r . , 7 .ft> , J ft , ft, ,ft , 1 -ft „ a .ft) «.*«(*) -ft) 

" /Oiu. <_ 4lv s) 4«- to ft* t ft />4*>- 1 4m. \) 4a*> to ft-* 1 ^ 

, /n , . 7 - ®) - "5 .ft) , 7 ft, .ft) (•> ft) . 7 -ft) . * fo) -ft) (O) ft) 

^ iOA~+- l ft— %) >iXw to Co -* c ft-*\) ft"*CO — c 4 vy- V ft-* c ft-? O ft— »io 

-ft) , •< ft) . .ft) ^ - ■».<•)_ *4 .ft). •» ftj> - 6 ). „ . . ft> -ft). . (•> 

*+ t »d^V )dwCO £#-*c — /QAv 7 t \) 4^ CO C*~?c 44^- • 

.7 ft> . (S) 'ft)”! S' - 7 -- ft) - , ft) » t ft) , ft), fo 

/fC- \) /T.— v Co ft* C J ’ ?4**> C 4—- S) 4-» CO ft-* 6 ft * W 

- *. 6 ) . 7 ft) , 6 >) -fe), ^,<t>, 30 „ - vftr- J 0 )*/' °-) 

— £ XU—- c /ft— \) /i. — CO ft-? C ft-* N> -f -*^ ft— C ft— \) ft— CO ft-* c ft-»N) 

. 2 s-.- ^ ^ ' ’/ /**, *fc>. ( 3 ) < - '*.,('?> . “? 6 > 

+ ' ^ ft— c ft— \) 4 m. CO ft—* C <fo* CO 4 m - L ft— \) ft-* c ft—? CO 

< 

5 


- J 5 *^ i ( '- , ^.v frf ? , 4 i-V W M* £<*^ «/>-*/*. w f *M? l ^ N)^ 


<,£>• 


4 T /,(•>- ,(o> ft . - r °) - /°' 

— ty C ft— «, ft-* t ftr» CO *■ C ft-* c ft-* CO 


ZpCc)* o 
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'• ??■**&. Tim* vr * w.^% . 


Element Q 

ft f" 


i. 5 ft 

"" ^ V C 


/■ \ ftr a) 't Gj -fij i- ' ^ .ft - (*j . (% w (•> *_ ." » - >-» •' 

^ ( Aj = e jh/^V l Wu - 3^Xlw V 4 wCo C«C O.V ^S-> 4 « -t- 3^ - 4 — V o*<- 

- ft . .6* . .ft 


.ft. 


(i * 

-r« 


-ft 


(oj -o , (o\ t«y £> £j 7 JL-. /•> .!•> .ft t ft -ft, 

t - 4 > ■^* 4 *<.U‘ Co->\) to ■> t o j ^ ' C< • 9 <- - 4 ««\J /<Cv Co Co^k. Co^»N) ' 


, j... . -ft ft . i , ft, -<«>, , ft 

4 ^ to*? t £.<v» v — ■ ,x— s "O Aw. to c-®<7 c C#-p Co 


/ 3 ftj - 


s' / \ C*>>f \f G-. , ,ft, ^ ,(P> £’ «A» .ft - , ft, /£ 

•> 5 l (_n.J = e A- ■jr’M 4 —c /W- to c*-* o +v cX*c 4 wU> ■+ -Tf 


//" - -ft . «) . o) 

' SX*S— L 4U.\) ^jT-r Co (. t*—?\) 


. -ft . (•) . 5 ft 


jfi> 


/.^, . -I*' . /•> . < AT» . .mV* • V . •■ Ml, 

— IOa**?~ t 4w>V /4-Cv. to Coo*. Co^sj — I /l-C*. V /CCw u> C<r-*c £o-»Cu 


.ft , a ft _ i & .ft 


/< . , ft , 9 A) ft , ft ft . . ft . ? .ft , ? ft - ift, . /» 

— - J - A * — a <!«-*» V* /lw CJ C*"? 4. Co-» V 4 /D4ft* - \? /4o>- CO C«r? C C*-> V 

+ /0 Co~ Coo C*J °'— ^ 4*^C^ > Co~ £**Urj CO *V cC.t Aw.'O A*ZmU) Llr>\ 

+&. H-*~Z ^U‘ 

Fur, Fc^uS*- [V.Vt CO^ I ’ 


cp;cj 


4 «o ° £*■* ^ 9 -f ^!»C 


- ft . ft ft. 


fpi 


/ O , - , w Ml 'Ml i- . — _ V-M , 

+ <^' 4 w.L C*^»t V co 4 cX~ l aC- V /SCrCo 


ft . *''**» , - , 

/Kv-Cj — c iC<‘ 7 V £-*^» ci 


.(oj 


^ /■ \ l£- m 'ft - '*> 

(j\) * - tf € 4 . /i— ^ to i 


, ft (®> 


. ft , /•> 


cfc^c 4ov. CO 



Element g? 

iC/coW'f 


2 «.Wlo/* C~J*J 


+ C ^'\f*AZ.cJ*c~,tJ*+ 2^.\>^J l l J q> 6 


sjcb' 0*^- 4 e (4j ^- , c f ‘- > 

£,k-) = 

-i*cC r, *i.v®Al.co 6 + t* t t' ft ^. v ®ar,w®- 3Xfc. , A^v®^ , u f:) 

-47 ^ V^\>®L-‘ v w ff ^ v®<>, to®-*- Jjhi*A®^\,«^*«U w fe 

- J^A 6^VVA®*uiAo®- , J* 

+ 36 •2^V\c_\ft^.' 1 cj^c*-, v®£», (J^-t-Sl^C r \u ^-J^AC. cJ’ 1 

- $L L^i <J»+ afirtfL, 
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- cJ -01 ^ v£4 J°< 

*-f J'> 
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Element M 


2 _ . - * ..ft) -2 - T.-ft> - 1 ft) ft) ft) 

C -* j «. -f ^5^ — t \J ,4-CC C-J £«-» to 


v /»_A .ftl 0 3 ft) 7 (e) 9 7 fci 

— £ /v — •y \J tT ■+ j -l*^- «T yj 

-% 0 ;^, 6 kVS 


VV , 1 ft). - , ft), ftj 

<**• V <)wU> 6 #-» 4 o" ■' 


,«) _ „ft) 
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Substituting Equation (D-3) into Equation (D-4) yields 


cd + 3Z 4iA.se® 
— r" v i*» 1 - 


(D-5) 


secular 


nonsecular 

The general resolution of v' w/ into secular and nonsecular parts can be indicated 

as 


,( 0 ) 


(5) (?) 

V = V, + v w 


(D-6) 


Equating Equations (D-5) and (D-6) results in 


V 


& 

* 



(4-23) 


Inspection of Equations (D-3) and (4-23) shows that and E^ have the same 
secular part. 


APPENDIX E - SOLUTION OF THE SET OF ORDINARY DIFFERENTIAL 
EQUATIONS HAVING t AS THE INDEPENDENT VARIABLE WHEN 
~ CONSIDERING EARTH OBLATENESS ONLY 


E. 1 INTRODUCTION 

The following set of ordinary differential equations to the order of e^°\ as 
derived in Paragraph 4. 3. 1. 1, will now be solved: 


i /•'_ 

d-t 


dco* = 
dt 


!s_ ^ A 

— ( T ‘ -•) 
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(? o) rt /„ «• . TlS A- .<*■>- © -.C 

a 1 /> + (f e ^‘ ' 1 ~‘“ 

(*-> .S) . fqt JL ,* - 6>// < - ^Y1 

— e c<k. l c co + c -4^ co (y~ c yj 
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I *i r o> (s> <9 ,<»> (Oj 

^ _ K> V /Zz Kz. /V» € C^-CO /, ^ „ 3.(5)’ 


dt: 


0-e®’) 




(4-36) 




d? 


d^/ e) 

d£ 


(l-e r °"y (l-e^y 


e®V - f> 


(5> 
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C* r «- ^ /, y . 2 . (o)\ 

) 3 A ( ~ ^ d ✓ 


(4-43) 


+ 




-y '4-*-*“ V Co ( / C ylj 
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The method of solution is based on one set forth in Reference 16. wherebv e (0 * 
( 0 ) 

and co are considered to vary simultaneously and the solutions for the remaining 
elements are obtained as functions of e (0 * and co ( °\ In consideration of this, the 
following (assumed) constants are defined: 
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Note that these constants are expressed in terms of the epoch values of the element 
functions, i. e. , ej^, and n^. The effect of this assumption on the 

accuracy of the resultant solutions can be minimized by periodically rectifying the 
orbit and updating the epoch values of the elements. (As discussed in Paragraph 4. 5, 
an "updating 4 rocedure" is used when numerically evaluating the solution equations. ) 

In terms of these constants, the differential equations become 
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Equations (E-10) and (E -11) will now be solved simultaneously. The solutions to 
Equations (E-12) through (E-14) will then be obtained as functions of e^ and 

E. 2 SIMULTANEOUS SOLUTION FOR e (0) AND w (0) 

To effect the simultaneous solution for e^ and the following transforma- 

tion parameters are used 
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Differentiating with respect to T yields (where the ’’dot" 

indicates 7 differentiation) 
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Substituting Equations (E -10)* (E-ll), (E-15), and (E-16) into Equations (E-17) and 
(E-18) yields after simplification 


i—c -n -I 


To the order of e, i. e. , ignoring terms of the order e (or smaller), these equations 
become merely 

<f ~ -*■ (E-19) 

^ •= (E-20) 

(NOTE: These equations agree in form with Equation (14) of Reference 16. ) 

In terms of the operator D = , Equations (E-19) and (E-20) can be expressed as 
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■ £ 

To solve these equations simultaneously (Reference 22, pp. 198-200), £ will first be 
eliminated. Substituting Equation (E-22) into Equation (E-21) yields 

c,c-l - cl •>, - y l (^) » o 

or 

(y'+cf)-') - c, c z 

As can be seen, this is a linear nonhomogeneous equation in It has the standard 
solution 

V = K, Cjt + « (E-23) 

where K and K are undetermined constants. 

X ^ 

To obtain the solution for £ , T) will now be eliminated from Equations (E-21) and 
(E-22). From Equation (E-21) 

* 1 - 

implying 
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Equating Equations (E-22) and (E-24) yields 
or 
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As can be seen, this is a linear homogeneous equation in £. It has the standard 
solution 

C 2 t + (E-25) 

where K and K are undetermined constants. 

”3 ~4 

Now, the constants K. , K , K„, and K have to be "adjusted” so as to make Equations 
(E-23) and (E-25) satisfy the original equations. This can be done by substituting 
Equations (E-23) and (E-25) into Equation (E-20) and seeking relations between these 
constants. Performing this substitution yields 

Ctr-> C-2^ 4- •+ Cl "t ■+• = O 

which implies 
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This equation is true only if 

= Kt_ 

and 

Kq = - K, 

Furthermore, it is convenient to express and as 
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where A (a positive number) and a are constants yet to be determined. Thus, 
Equations (E-23) and (E-25) become 
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Since, from Equations (E-15) and (E-16) 
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it follows that the desired solutions are 

e'®= [f \ 2 A *:~(c*V+4 + ( %)* ] 1/1 (4-33) 
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It now remains to determine the constants A and a. This can be done by evaluating 
Equations (E-26) and (E-27) at the epoch time_t resulting in 

£ 0 * A t+~-> C.(-x (E-26) 

°7 0 - Pi 4 O + (E-29) 


Multiplying the first equation by sin (C t + a), the second by -cos (C t + a), and 
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then adding yields 
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which becomes by Equations (E-15) and (E-16) evaluated at 
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Then Equation (E-30) can be written as 
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Multiplying Equation (E-28) by cos (C^ + a), Equation (E-29) by sin (C,^ + a), and 
then adding yields 
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E.3 SOLUTION FOE i (0) 

By Equations (E-15) and (E-26), Equation (E-12) becomes 
-^4r — dy A c*~7 C(- -l ?■ -#■ o<3 

which yields upon integration 
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where K is the integration constant. Substituting Equation (E-27) yields 
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The constant K can be determined by evaluating this equation at the epoch time^ 
resulting in 
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Hence, the desired solution is 
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E.4 SOLUTION FOR O 
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By Equations (E-16) and (E-27), Equation (F.-13) becomes 
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which yields upon integration 
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where K is the integration constant. Substituting Equation (L-26) yields 
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The constant K can be determined by evaluating this equation at the epoch time^, 
resulting in 

*= cof 

Hence, the desired solution is 
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E. 5 SOLUTION FOR M (0) 
From Equation (E-14) 
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This equation can be linearized by approximating e^ by its epoch value ej^ (referred 
to as "backlining" e*%, resulting in 
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Then, from Equation (E-16) 
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where K is the integration constant. 

From Equation (E-27) 
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Evaluating K at epoch time^ yields 
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As mentioned in Paragraph 4. 3. 1. 1, the Keplerian change in M taking place during 
the time interval (t - t Q ) must be added to the above equation. 
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APPENDIX F - COMPOSITE SOLUTION OF THE SET OF ORDINARY 
DIFFERENTIAL EQUATIONS HAVING T AS THE INDEPENDENT VARIABLE 
WHEN CONSIDERING EARTH OBLATE NESS AND ATMOSPHERIC DRAG 


F. 1 INTRODUCTION 


The following set of ordinary differential equations to the order of e^, 
representing the combined effects of earth oblateness and atmospheric drag 
(Paragraph 4. 3. 1. 2), will now be solved: 
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Again, the solution procedure will be to solve simultaneously for e^ and<t/°\ 
employing the same transformations and similar approximations as in Appendix E. 
These solutions will then be utilized to obtain solutions for the remaining elements. 
In consideration of this, the following constants (in addition to those given in 
Appendix E, i.e. , C - CJ are defined: 
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In terms of these constants and those defined in Appendix E, the differential 
equations become 
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Note that these equations are extensions of the oblateness -only equations (Appendix E) 
and are obtained by adding the drag perturbations to the respective oblateness 
perturbations. 

F. 2 SIMULTANEOUS SOLUTION FOR e (0) ANDgX 0) 


Applying the transformation equations 
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to Equations (4-71) and (4-73) results in 
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These equations can be linearized by approximating e' ’ by its epoch value e 
(referred to as ,r backltning" e^) f resulting in 
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Differentiating the first equation, substituting from the second, and employing 
the differential operator notation of Appendix E yields 
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The characteristic solution for this nonhomogeneous second-order equation is 

given as „ 

^ b ^ r ^ * 1 

£ c = ^,ux>ct + c t J 

where and X are the integration constants, and the particular solution is 

e ^JSL 
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Hence, the complete solution to Equation (F-21) is given as 
£“ e *e + 

The solution for is similarly obtained as 

'’1* (F-23) 

The integration constants X and X can be obtained by evaluating Equations (F-22) 
and (F-23) at the epoch time resulting in 

*>“ e v"C(£°+ ^77*) -<S. + (F-24) 

+ (F-25) 

From Equations (E-15) and (E-16), the desired solutions are 

(£*+• (4-72) 

W® = t2~*' <4-74) 

where and J7_ are given by Equations (F-22) and (F-23), respectively. 

Recall that Equations (F-22) and (F-23), even though exact simultaneous solu- 
tions to Equations (F-19) and (F-20), are still approximate solutions to Equations 
(F-17) and (F-18) since these equations were linearized by backlining e^\ The 
validity of this approximation is demonstrated by the data in Table F-l, which con- 
sists of time-point comparisons between numerical integrations of Equations (F-17) 
and (F-18) and numerical evaluations of Equations (F-22) and (F-23). The data 
span a 60-day time period with an integration step size of 128 seconds and a print/ 
evaluation step size of 0. 5 day. The orbit eccentricity is 0. 036. 

In Reference 16, plots of the motion in the £, y -plane are given when considering 
oblateness only. The rj_ solutions presented there are of the same form as those 
derived in Appendix E, i. e. , 
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Analytical 
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Analytical 
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-0. 1707371 7xl0 _1 
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0. 13638160xl0 _1 

10 


-0.32363856xl0 _1 

-0. 32367086x10" 1 


-0. 9799636 7xl0~ 2 

-0. 98002792xl0~ 2 

20 


0.41330305xl0" 2 

0. 41331292xl0~ 2 


-0.31283737X10" 1 

-0. 31290135x10" 1 

SO 
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40 
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50 
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0. 1344061 7x10 1 

60 
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-0. 18930974xl0 -1 

-0. 18943342x10" 1 


Table F-l. Comparison of Numerically Vs. Analytically Integrated Values of v 
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Plots of these solutions are circles centered at (0, C/K)„ The magnitude of 

e^ (t) is given by the length of the line segment from the origin to the point 

(£(t), rj(t)) on the circle, and the angle turned by this line (relative to the + £ 

(0) 

axis) measures oj (t). 

Since drag perturbations tend to diminish e, it was expected that the correspond- 
ing plots for £ vs. Tl, as given by Equations (F-22) and F-23), would depict this 
decaying effect while retaining the basic characteristics of the oblateness-only 
plots. Figure F-l reveals that this is, indeed, true. The plot shows a 360-day 
variation of £ with for the orbit of 0,0364 eccentricity. 


The solutions for£, T]_ given by Equations (F-22) and (F-23) represent the 
combined effects of oblateness and drag. Although these solutions are functionally 
different from the corresponding solutions for oblateness-only (as given by Equations 
(E-26) and (E-27)), they are identical in the limiting sense, i.e., as the drag 
perturbations tend to vanish. This can be seen by allowing b (or D ) to vanish in 
Equations (F-22) and (F-23). 
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*1 ANALYTIC SOLUTION 


Figure F-l. Plot of £ vs t? as Given by Equations {F-22) and (F-23) 
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which reduces to the oblateness solution 

F.3 SOLUTION FOR i <0) , Q (0> 

Since drag does not affect i_ and 0_ (at least when the drag model is tangential), 
the solutions for these elements remain the same as in the case for oblateness- 
only and are given by Equations (4-37) and (4-39), respectively. 

F. 4 SOLUTION FOR B (0) 

The solution for is obtained by approximating e^and B^ by their epoch 
values on the right-hand side of Equation (F-15). The resulting equation can be 
directly integrated to yield 

r*= (4-80) 

(NOTE: Numerical results have verified the accuracy of using the backline epoch 
value to effect the integration. ) 

F. 5 SOLUTION FOR M <0) 

Equation (F-16) can be rewritten (using Equation (E-15)) to yield 

= C( " ~ ^ c i~ **) y l + <£ 0 + i| e^J 

The epoch value of e^ , i.e. will again be employed, and the solution 
form is given by 

where K is the integration constant. 
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From Equation (F-23) 
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bt e „ 

“ b T7^7 {( ?. -4 (S r l? 4 7\,c) 2^7^ 


( 4 - 83 ) 


where the integration constant K has been evaluated at epoch time t^. As mentioned 
in Paragraph 4. 3 . 1 . 2, the Keplerian change in M_ that takes place in the time 
interval (t-t^> must be added to the above equation. 
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APPENDIX G - PROCEDURE FOR DETERMING FUNCTIONAL FORM OF C(t) 


G. 1 INTRODUCTION 

( 2 ) 

An outline of the procedures necessary for obtaining the E solutions was 
given in Paragraph 4.3.3. Inherent within these procedures are the steps necessary 
for determining expressions for C(t). In this appendix, the details that must be 
considered in implementing these steps are discussed. Each step is illustrated 
by working out the solution of C(T) for the inclination^, denoted by C.(t). Because 
of the complexity encountered in solving for C(t), a number of simplifying assump- 
tions are made. However, the final expression for C.(t) is shown to agree re- 
markably well with computed values from the GENPUR program, thereby indicat- 
ing that highly refined equations may not be necessary in determining the functional 
form of C(t). 

G. 2 SECOND-ORDER EXPANSION OF THE BASIC DIFFERENTIAL EQUATIONS 

2 

The first step in the procedure is to expand, to order e , the basic differential 
equations of perturbed motion. These equations for the gravity perturbation 
are Equations (2-34) through (2-39). Consider, for instance, the portion of the 
equation for i_, Equation (2-36). 

jiL = - ^/^ e 1 fi + e c„ V ; ^.3. (2 ' 36a) 

The required expansions may be obtained using the methods of Appendix B, i.e. , 

6 7 = 6 7e fl4 eP 

(it - (t -+ 6 feO 4 

0-e.')- 7/2 = (I - ^ (G-l) 

/J X*. 0.L ~ 2 C® ‘ +■ 6 2 Qi ^ 
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These expanded expressions can be written in the form of binomials 

V'= X, + <EVi 

(G . 2) 

(t-e')' lb = Xs+fi-f, | 

C2t — )(i^4 6 v |i| 

-t 

where 

X, * 6* 1 " 7 

x a = 

X^- = 

and y, - 7 6® ‘gffl 

Vi = 

Y, - 7 ((-«**)-* e ft «® 

Y ¥ = 

Vs- - 2J i c~,2j°> 

The basic differential equation expanded to order €, in terms of these binomials, 
is merely 


= - £ 6 /XT r t Y*. & *5 X, X f/ ) 


(2 -3 6b) 
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(G-3) 


The same equation expanded to order £ is 

* *“ "2" ^ 1 fa Xi X 3 X<^- ^ £-V>^ <c* (y, X-j XjXy Xj'+ X, V2 X ? Xy 

+ x, x 7 v, x v x y 4 x, x 7 x 5 Vy + x,x,x,x„y*.) 

or, in more concise notation 


"dr - £ £ + 


which is equivalent to Equation (4-115). 


tic 


(NOTE: Equation (G-3) illustrates a convenient method for expanding the basic 

2 

differential equation to order £ . First, convert the expansion of each term 

to the binomial form of Equations (G-2). Then, insert the y part of each binomial 

2 

into the differential equation, one-at- +-time, to get the e expansion.) 

The complex nature of Equation (G-3) is veiled by the notation. For instance, jr is 

Vi - 76'f K/e)fj:/e)- 4 c/F)] 

+ ^’t- fc ^’v rc ’- + £ e ("£)\ 


The remaining ^ terms are equally complex. Some simplification is possible, 
depending upon accuracy required, by neglecting terms containing certain powers 
of e^ and by holding certain elements constant, such as i_. (Use of the FORMAC 
program would allow retention of these terms in the expansion. ) 

One assumption that will be made at this point in working cut the solution for 

the element i is that all terms involving powers of e^ will be ignored (i. e. , e^, 
(0) 2 

e v ’ , etc.). Thus, 

- ( /-e ) = l.o 

N*' 0 
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The L and functions become 
£*£*■- "5 ' *"c ^X# X? X v Xy) 

^i--j[^M ^(i, X*X V %g + X, V 2 -#• X, x 2 y s + x, x t x v y s ) 

where, recall, f is the expansion of the basic differential equation to order e, 

■* 2 

and g. is the expansion to order £ . To further simplify for later analysis, let 

d,= y, X 2 X V X^ 
s V** ^ 

d H = X,X 2 V*<X^ 

* X/X^XyV^ 

so that 

r”c X (clj + d a -f -f 

G. 3 PARTIAL DERIVATIVE OF E (1) 

The next step in the solution process is to compute the partial derivatives of 
the E^ solutions with respect to £. Recall from Equation (4-28) that these solutions 
have the form 

H^C^ey ijeU‘>i ]+...+ c e m 

or merely 

B r 'K C e (t) 

The partial derivative is then 

41 ?- ^ 

HT~ 

so that the crux of this step is in determining dtJdT. For those solutions 

N 

that have no secular terms in I (E) (or I (E)), such as _B, this step is relatively 

u u 
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straightforward, but tedious. The element i_has no secular term in I (i). 
Neglecting T terms, its super-one solution is 

O 

K t (i)i,(0 * ojt) 


such that 

(„(:)= K-iC 


The partial derivative of f N (i) is merely 


a ? 


a Cc) CO] 


One question that arises at this point is whether 


iO- e® 1 ) 1 

can be considered constant. The term is constant if gravity perturbations 
only are being considered, and can be treated as constant for fairly long intervals 
even when drag is present. The e^ term, however, has an important long- 

periodicity due to It has not yet been determined whether the e^ variation 

- ^ ( 0 ) 
with^ should be included. However, since the example is neglecting e terms 

and is considering only gravity perturbations, K (i) will be treated as a constant. 

Evaluation of 31 (i)/3t requires further analysis. I (i) is composed of 
terms such as 

“* #4^ ml L CO CO V 


so that the partial derivative involves determining de^°Vdt, di^/dt, do/^/dt, and 

d J°Vdt. The first three of these derivatives are available from Equations (4-32), 

(0) . — 

(4-36), and (4-34), respectively. Determining dv /dt is more complex, however, 
since the explicit expression for is not immediately available. To obtain such 
an expression, the Fourier-Bessel expansion may be used, i.e. , 
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V^ o) = Ao^+ 


& • • 


The ds/°V dT derivative can then be expressed in terms of dM*°Vd?and de^/dt' 


dv r#J 




(G-4) 


If only the X gravity perturbation is being considered (de*°Vd? = 0) and terms 
, (0) 

involving e' are neglected, then 

d V* ' _ ^ *■■><•) 
a? ~ 1F~ 

An expression of dM* 0) /dt is provided by Equation (4-43). Notice that Equation (4-43) 
does not include the Keplerian variation of M*°\ which is dependent on 7 rather than t. 

Even with these simplifications, the partial derivative of only one term of 
I 2 (i) is very complex, i.e., 


~ ~ a 1 ( ° )/tZ ^ 12 v fo) ~ 2 2 2 \> (o) 

- e (o) ^2: Co) ^ Qc 

Substituting Equations (4-32), (4-36), (4-34), and (G^) yields 

- a e r<I L, 

- ^ / )2 rYa- 4 

0- ?')-'(% cJ'^L f U K 

(I- _ e<V. 2e®o, 

( i- e^T’V/- *^ ft ) - ^y sh - 

£3 c'V. J->(i- <J->(i- *aL.'l toyjJ 


deft 


+ derivatives of 8 more terms. 
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As noted, this is only the expansion of one term in nine. Thus, evaluating the 
derivative of even the simple type 1(E) , which has no secular terms, quickly 
becomes very involved. 

To return to the more general form of an I (E) integral, which contains secular 

( 1 ) 2 

terms also, consider the equation for Q , i. e. , 

jx <r> = KaMfcOd- K ] + . . . ♦ C^Ci) 

The terms within brackets are an expression of nonsecular (with respect to_t) 
parts of an integration. For the element j_, there were no secular terms in 
I 2 (i) so that S 2 (i) was equal to zerc. For the element Q, however, the I 2 (0) 
integral contains a secular term, - 1/2 \f^ cos i^. To remove this secularity 
from the brackets, S (O) is 

SaM* 

Assume for the moment that 1(0) was composed only of the secular term, and call 
itl'(fl). Then 

To compute 3f N (0 )/8t, the partial derivative of the right-hand side of this equation 
must be evaluated. Substituting Equation (4-45) for yields 

(S'-J'i) = . . .-,/*>* (G-5) 

Next, using a functional form of Equation (4-44), i.e., 

a/ o) = a + (4 -44a) 
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Substituting Equation (4-44a) into Equation (G-5) 

Cv^-a - c "'t)= 2n 6 +..,-A 

- /*> r, f© + 2e r Vv.r^ ft * 4 2 (G ~ 6) 

and taking the partial derivative of Equation (G-6) with respect to t_ 

■■jb (iP-~Pt) = 0 * 2e a a,/->®) -+^2^-^^+... 

at - dc 

which is the same as Equation (0-4). Thus, the presence of secular terms within 
an I (E) integral presents no new problems, but merely adds similar terms. 

Performing the differentiation of all terms of each E^ solution with respect 
to f and making the necessary substitutions will result in extremely long expres- 
sions. (The need for an automatic manipulation language such as FORMAC in this 
step is indeed evident.) 

To continue with the example for 2 , neglect terms of Equation (4-89) of order 
(O' 

e ' (and higher) and omit the T terms 

+ 2 ^ + C L (£) 


Assume also that 



O 


dKj^O 

d* 


- O 


The3f N (i)/9t derivative is then 
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= K x 2^3;V^\ Ca) A^2u> r °>] 

<G-7) 

-( 4 ^ 2v ro) + Qcj>~] ^“’j. 

(To facilitate work in the next step, expressions for do/°Vdrand dJ 0) /d?will not 
be inserted at this point. ) 

G.4 SECULAR PARTS OFT INTEGRALS 

A second application of the first uniformity condition requires determination of 
secular parts of the integration of g and af^/at with respect toT. However, g and 
8f N /8t are primarily functions of v (0) rather thanT. Thus, It Is more convenient 
to change the variable of integration from! to by the relation 


T* e (° >co^> v <**) z 


The following integrals must now be evaluated and the secular parts determined 


fv*-k 


1 J 1 ( V 


/ 


l£ <lc = J SF 


iC: 6<°>Yi-, dv 


This step is straightforward, but very tedious if no simplifying assumptions are 
made. FORMAC has been modified to evaluate these types of Integrals and to 
then extract the secular terms. 
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To continue with the example for the element^, consider the integral of 

+ d^-4- d/) 

There are some terms common to the which may be extracted, i.e. , 

% ~ <f 7 ' 6 t<5> C ( he (0> c^ \> <Pi y(? l + -i H: s ) 

where 

t,* (<> 6 \ C^v^- e^v^v'’") 

Z 2 = 7 Q c *(h t r "V,v '*> ) ^ a ‘‘•Ux. 2c *> 

z 4 = 2 €><\- <Y/ -« e v/’ 1 ) ^ 2 ^ 

2s . = 2 e C °'t*-> v <S ')' 4 “- 2 ; c °> u,Jk S J 

Then, 

fct d* = - * Is) df (G-9) 

(0)2 3/2 

Neglecting the (1-e ) term in Equation (G-8) and substituting into Equation (G-9) 

yields 

= - -£ f * 6‘° |3 jVz,+ Zj-* ? h ■*■ if] dv^ 


To evaluate integrals of each z. term, without the use of FORMAC, all terms 

I 

containing e^ and higher powers of e^ are ignored. The expression for z^ 
becomes 


2,= 36 %£~2; f, iiz.3S , e n cr-v fv 
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In obtaining the value of e^' to be used In z^, only the J portion Is considered, 1. e. , 

- io'V, Jjf-AiJt 

-4*-’;^’^^- c^,v ro) J 

+ C e £«■) 


The expansion of cos sin 2u^ is 


3*4* - 2 c^t 0c4 i /iL*'J' i <~4'J'‘ > +/C**. 0 cJ“ i c*-> 4‘L 2 V*V'. 1 


Hence, 


i 

f2^®AC.v r -k-?^ ^-^c'^k-’u^co^vf"- 
2j’ J - v®- 

+ ^-l 
B»V J 


Hereafter, only those terms which are secular In v , after Integrating with 
respect to , will be retained. Such terms are of the form cos^v^ sin^v^, 
wbere M_ and N_ are even integers (or zero). Thus, the secular terms in z are 
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t, iec ^ ro f<^£«2u ll, *. 


2 JV-2 J , Wv*V-V*'-,V , V 2<-.'-^. V '<- V-' 


2 . t.(b . o ? (•) 


<- 4 « » . 


w 


) . (o) 3 2P - i .£) .i (•) - fo> . ij (•) 2 (i > 

■+ >W 2u> J C*-+ V — > c A*** to /iA^w'u ) Ctr 7 \j ' 


4- .a.4A.V^2c- t,> ^V i k- , u W - 2^?~‘Wv |6, , 4 i.V*'J 


Integration and simpllcation yields 

fs, df =• -f- r. 0<J‘XQ- ?*_’<: s, Jv ^ 

' sec 

By similar procedures, it is found that 

K« i€ = 0 
K« «** - ® 

where in the ^integral u^ = v (1) + a/ 1 *. The secular part of the 7 integration ofj^ 
is then 

fat tit * - r c M 6*^ 2: C ^'2^( 2- 3a^C (G-10) 

itt 

Next, consider the integral of Bfj^/ai with respect to 7. Assume in Equation 
(G-7) that dcJ°Vdt and d J^/dT are constant with respect to a t_ integration . 

Referring to Equation (4-34), it is obvious that the secular rate of change of 
with respect to Tdoes not depend upon T. It is not so obvious for dv^/dtl 
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since v 


( 0 ) 


has a Keplerian variation that depends upon”t . 


In functional form 


+ f (Z) 


However, its derivative with respect to_t is not a function of t_, i. e. , 


d l'O ) 

"J 1 \ 

& - + (t) 


Consquently, ds/^/dt is constant for the t_ integration. 

It is again convenient to change the variable of integration from T to by 
using Equation (G-8), where the e^ term is neglected, i.e. , 


<{ t ~ 




<i^ 0) 

e (o) * 


The integral is then 

f I X _ Kl (i) {<4*, /a. , /J.b) . ^ (el . 'i f») . - _ to) 

J 0(o)T J j A^r. v — 2 Ur-> 2 


(o) ~ fco) ioj 

GO A \) C^9^J { J 


)<w° 


<4? 




Eliminating all but the secular terms gives 


ra^u) = 


J j.? 

SEC 


w 


d? 


(G-ll) 


G. 5 INTEGRATION WITH RESPECT TO t 


The differential equation for C(t) from Paragraph 4.3.3 is Equation (4-122) 


< K;fa r 

dt 


- f v & - 




(4-122) 
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The final step in the solution procedure is to solve this equation for C(t). Sub- 
stituting Equations (G-10) and (G-ll) into Equation (4-122) and then integrating yields 


( (O . Go . o . o> . o. ») d<Y> ° t-r 


To evaluate the integrals, substitute 


N . (° __ (“i r 

V — /v-t XT 


«,ro- 


so that 




- (i 6" 




(G-12) 


At this point, the E^ solutions for each element (Equations (4-33), (4-35), (4-37), 
and (4-44)) should be substituted into Equation (G-12). The resulting integrals will 
be complex and very difficult to evaluate. For instance, the sin 2 co^ term would 
have the form 


/ 2 1 ,~- 1 1 

P\ C^-~3 (C-?^L •+‘*0 ^ J 


which must be combined with other functions of_t and then integrated. 

For the inclination example, a different approach was taken. All elements 
within the integrals were assumed constant, except to Then, from Equation (G-12) 

£<-(?>* - £ 6“'^. 3;*Y2- W<- '* V*'p^ de 

2co'”’ jV* 1 
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In the first integral it is easier to express dT in terms of d <J^ than to express 
sin 2 <J^ in terms of t. Using the portion of Equation (4-34), 

Ct 




dco (0 


J°> (2- *>) 


In the second integral, the dts cancel. Thus, 

Cc (t)-% (U. 2 J'> Aj* 


(2- gf *~*Z*') 


Integrating and collecting terms gives 


CiCt)* -k ) + l"] «f< (G-13) 


where K_ is a true constant of integration. Equation (G-13) is the final expression 

for C.(t). 
i 

To test the validity of this solution, it will be compared with results from the 
GEN^UR program. Appendix H presents solution components and associated 
constants as computed by GENPUR. 

The plot of eC^t) is reproduced from Appendix H and shown in Figure G-l for 
the low-eccentricity orbit case. Evaluating the constant K_ and inserting proper values 
for the initial conditions in Equation (G-13) results in 

CiC *0- 0.27 o-7 2 <J* (>*<0 


In Figure G-l, the values of C.(t) have been multiplied by £and expressed in degrees. 
The corresponding solution equation would be 

o.osg’ 
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Figure G-l. CC.(t) for the Low-eccentricity Orbit 



The dashed line in Figure G-l is a plot of this equation. It matches the GENPUR 
solution extremely well. 

In summary, results of the simplified analysis to obtain the functional form of 
C.(t) are very good. These results bring up the question of how much improvement 
is really needed or desired, especially since a great deal more effort would be 
required to eliminate the simplifying assumptions. 
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APPENDIX H - PLOTS OF SOLUTION COMPONENTS AND 
ASSOCIATED PARAMETERS 


The second approximation solutions to the equations of motion result in ex- 
pressions of the form 


E (t) = E (0) + € E (1) 


where E is any element in the set (e, i, oo, B, or M). The E^ solutions are mean 
elements and contain the secular and long-periodic variations. The E(t) solutions 
(i.e. , E^ + € E^) are osculating elements and contain short-periodic in addition to 
secular and long-periodic variations. This appendix presents plots of these solution 
components for each element when considering a low- eccentricity orbit. The initial 
conditions for this orbit are shown in Table H-l. 


Asymptotic series solutions were first derived when considering only the oblate- 
ness perturbation (J^ and first-order effects). Later, combined solutions for 
oblateness and drag were obtained. Therefore, two sets of plots are presented in 
this appendix; one for the oblateness-only solutions and one for the combined solu- 
tions. In both sets, Brouwer's equations are used for second-order secular effects 
of «J and J, on oo, and M. 

Deeply involved in the solution procedure are assumptions that various param- 
eters remain constant, at least relatively so for short intervals of time (2 or 3 days). 
It is known that some of these parameters (the C(t)'s) contain important secular and 
long-periodic variations. Other parameters (C. through C ) remain fairly uniform, 
at least for the oblateness-only solutions. Plots of these parameters are also pre- 
sented in this appendix. 

Figures H-l throu Ji H-12 show the solution components and associated param- 
eters for the oblateness-only solutions to the low-eccentricity case. Figures H-13 
through H-24 show the corresponding solution components and associated parameters 
for the combined oblateness and drag solutions. Interesting differences can be seen 
in the behavior of some of the parameters when drag is added to the solution. Many 
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Mean Elements 


B 

-1/2, 
(km ) 

0.012340 14 

a 

(kn) 

6566.5731 

e 


0.00559414 

i 

(deg) 

50.01120 

n 

(deg) 

152.47131 

Ui 

(deg) 

52.62626 

M 

(deg) 

2.91685 


State Vector 


X 

(km) 

-4872.6530 

y 

(km) 

-1364.7471 

z 

(km) 

4124.7041 

X 

(km/sec) 

4.41679 

y 

(km/sec) 

-5.51029 

z 

(km/sec) 

3.39451 


Pate 

April 1 , 1971 
C D (A/m) 

2 * 

0.0002 m /kg 

Table H-l. Initial Conditions for the Low -Eccentricity Orbit 
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of the parameters are functions of B^, which is constant for oblateness only but has 
a secular variation when drag is added. Thus, the parameters take on a secular 
variation which was not present in the oblateneps-only solutions. 

In the plots of the solution components, the secular and long-periodic trends 
are easily recognized. However, since these are plots of points at 6-hour intervals, 
short-periodic trends cannot be distinguished because they appear as somewhat random 
fluctuations about the mean. 

The units of the element solution components depicted in Figures H-] through 
H-7 and H-13 through H-19 are as follows: 

e - unitless 

to, i, £2, M, V -deg 


The units of the constants and functions depicted in Figures H-8 through H-12 
and H-20 through H-24 are as follows: 

C r C 2 , . . . , C ? - rad/hr 

cl - rad 

C (t) - unitless 
e 

cyt), C.(t), C 0 (t), C M (t> - deg 
C 3 (t> - km' 1/2 
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Figure H-l. Eccentricity Solution Components 
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Figure H-2. Argument of Perigee Solution Components 
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Figure H-3. Inclination Solution Components 
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Figure H-4. Ascending Node Solution Components 
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Figure H-5. Reciprocal Square Root of Semimajor Axis Solution Components 
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Figure H-6. Mean Anomaly Solution Components 
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Figure H-7. True Anomaly Solution Components 
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Figure H-8. C , C , and C Constants 
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Figure H-ll. f C^(t), €C (t), and e C^(t) Functions 
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Figure H-12. € C fi (t) and c C^.ft) Functions 
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Figure H-13. Eccentricity Solution Components 
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Figure H-14. Argument of Perigee Solution Components 
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Figure H-15. Inclination Solution Components 
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Figure H-16. Ascending Node Solution Components 
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Figure H-17. Reciprocal Square Root of Semimajor Axis Solution Components 
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APPENDIX I - DISCUSSION OF OSCULATING AND MEAN ORBITAL ELEMENTS 


There are two basic types of orbit elements , osculating and mean. A set of 
osculating elements represents an exact definition of one point on the orbit and, there- 
fore, is equivalent to cartesian coordinates of position and velocity at that point. 

Though osculating elements may often be difficult to obtain, there is no ambiguity in 
their definition. 

Mean elements, on the other hand, are not so well defined. As a satellite 
travels in its orbit around the earth, the osculating elements at each point in the orbit, 
when plotted versus time, will exhibit periodic fluctuations. These fluctuations may 
be either long-period (of the same period as w) or short-period (of the same period as 
J/). A mean element, as the name implies, is an expression for the average value of 
the osculating element (without periodic fluctuations). Usually, a mean element is 
defined as the osculating minus only short-period fluctuations. 

In engineering activities other than astroi^mamics, the standard working elements 
are osculating elements; for instance, the powered flight trajectory analyst will usually 
present insertion conditions in terms of a position vector and a velocity vector. 

Various theories, therefore, have been devised to compute mean elements from a set 
of osculating elements. Some of these theories have peculiarities in that their defini- 
tions of mean elements are tailored toward use in a specific general perturbation 
theory. Consequently, these elements are not truly mean elements but, rather, 
starting conditions for that particular theory. The definition of ’’mean’’ semimajor axis 
by Kozai (Reference 23) is a good example of this type of definition. Assume that an 
initial osculating value of semimajor axis, a, has been provided. To compute the mean 
semimajor axis, a, according to Kozai, first subtract short-period fluctuations: 

<x 0 - a - cja^ 


where 


cU* - Ja {%( I- l e< ) h ] + (r)^ i C+* 2 (v+0o ) ]■ 
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At this point, one would have the standard value of mean semimajor axis (oscu- 
lating minus short-periodic). However, Kozai continues and "conveniently" defines the 
mean value as 

5= a .ft- <i-i> 

which he needs for use in his general perturbation theory. The equations presently 
being used within the MSFC Orbit Lifetime Program (Reference 20) require the Kozai 
a element. Moreover, it is believed that the "mean" elements given by SAO in their 
reports on past satellite histories (Reference 19) involve the Kozai a definition. The 
GENPUR program, on the other hand, requires use of the a- type definition, as do 

V 

most other general perturbation theories. It might also be noted that output of the 
MSFC transformation program for "mean" elements uses a Q as the definition of semi- 
major axis. 

Another theory which has been used to compute mean elements is that due to H. 
Small (Reference 24). For the sake of identification, results of his theory have been 
termed "smoothed" elements. Small's theory involves removing short-period fluctua- 
tions from the radius and velocity vectors. If r, v , and v are Initial va«ues of the 

XV Li 

osculating radius . cadial velocity component, and normal velocity component, respec- 
tively, then the mean values are 

r = r - cic 

These mean values are then used in a standard coordinate transformation pro- 
cedure to compute mean (or rather "smoothed") elements. For instance. 


where 


M * - V + v t 2 

In Reference 25, it was found that the essential difference between the Kozai a^ 
and the Small a is 



The "smoothed" elements of Small are those used by the MSFC Orbit Lifetime 
Program for integration. However, it must be remembered that the A-E-P (Reference 
26) equations within the Lifetime Program require the Kozai ji given by Equation (1-1). 

Other methods for computing mean elements are given in Reference 26. Initial 
mean elements for Brouwer's equations are computed by the following procedure. 

Sets of osculating elements over some time interval are required. The long-periodic 
and short- per iodic variations are computed for each element in every osculating set by 
using Brouwer's equations. Then the secular motion is computed referencing every 
set to an arbitrarily chosen epoch time. Mean value sets are obtained by subtracting 
the secular and periodic terms from each respective osculating element set. These 
sets are then averaged, yielding one initial mean element set. Notice that, if only one 
set of osculating elements is available, the procedure is similar to that of Kozai. For 
instance, initial mean value of semimajor axis is 

a = c, - exp 


where 


<X P = fz <%<=’)“ 5/2 ] 

-+ " 5 ( I- - + > 2 


The Brouwer expression for a^, in fact, is equal to the Kozai expression for 

da . However, there are no further steps in the Brouwer procedure, so his mean is 

s 

truly osculating minus short-periodic terms. 
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The importance of having well-defined mean values cannot be overemphasized. 
For example, it has been found that an initial error of only 0.4 km in the initial mean 
value of semi major axis can cause a 360° error in mean anomaly at the end of 330 
days. 
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